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FOR  OCEAN  ENGINEERING  APPLICATIONS 


Abstract 


This  report  discusses  the  analysis  and  identification 
of  nonlinear  system  dynamic  properties  by  stochastic 
techniques,  on  measured  experimental  data.  Two  ocean 
engineering  applications  of  this  material  are  developed 
9f  CQncern  to  NCEL  representing:  (a)  nonlinear  wave 
force  problems,  and  (b)  nonlinear  drift  force  problems. 
General  models  are  formulated  consisting  of  parallel 
linear  and  nonlinear  systems  where  the  input  data  can 
b.e  Gaussian  or  non^Gaussian.  Formulas  are  stated  for 
Statistical  errors  in  estimates  from  measured  random 
data  to  help  design  experiments  and  to  evaluate  results. 
The  last  section  of  this  report  derives  various  useful 
Input/ output  relationships  when  stationary  random  data 
pass  through  three  types  of  nonlinear  nonsymmetri ca 1 
systems  of  physical  interest,  I 
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NONLINEAR  SYSTEM  STOCHASTIC  TECHNIQUES 
FOR  OCEAN  ENGINEERING  APPLICATIONS 


INTRODUCTION 

The  Naval  Civil  Engineering  Laboratory  (NCEL)  has  sponsored  develop¬ 
ment  of  a  new  approach  for  random  data  analysis  and  identification  of 
nonlinear  systems  under  NCEL  Contract  N62474-86-7275  with  0.  S.  Bendat 
Company,  the  contractor.  This  research  has  been  conducted  in  a  series 
of  tasks,  each  of  which  resulted  in  an  informal  report.  During  this 
period,  work  on  related  matters  was  performed  by  the  J.  S.  Bendat 
Company  for  other  sponsors,  and  much  work  was  also  done  as  independent 
research  and  development.  The  basic  mathematical  aspects  of  these  new 
results  is  now  being  published  in  book  form  by  J.  S.  Bendat  with  John 
Wiley  &  Sons,  New  York.  However,  some  of  the  material  that  relates 
specifically  to  ocean  engineering  will  not  be  included.  This  present 
NCEL  contractor  report  complements  the  book  by  emphasizing  certain 
matters  that  represent  important  new  practical  ways  to  analyze  and 
identify  nonlinear  system  dynamic  properties  from  measured  data,  and 
by  focusing  on  ocean  engineering  problems  of  concern  to  NCEL. 

This  report  discusses  the  following  topics: 

Section  1.  Parallel  Linear  and  Nonlinear  Systems. 

Section  2.  Nonlinear  Wave  Force  Models. 

Section  3.  Nonlinear  Drift  Force  Models. 

Section  4.  Nonlinear  Nonsymmetrical  Systems. 


Techniques  are  explained  in  Section  1  on  how  to  replace  a  large 
class  of  single-input/single-output  nonlinear  models  by  equivalent 
multiple-input/single-output  linear  models.  This  is  a  significant 
achievement  because  it  permits  complicated  nonlinear  problems  to  be 
solved  by  known  linear  procedures.  Input  data  to  the  nonlinear 
systems  can  be  Gaussian  or  non-Gaussian;  the  data  can  be  stationary 
random  data  or  transient  random  data.  Formulas  for  non-Gaussian  input 
data  require  the  computation  of  conditioned  spectral  density  func¬ 
tions.  Formulas  for  Gaussian  input  data  require  only  ordinary  auto- 
spectral  and  cross-spectral  density  functions  when  the  parallel 
nonlinear  systems  are  square-law  and  cubic  systems.  Square-law  and 
cubic  nonlinear  systems  are  the  nonlinear  systems  that  occur  when  one 
obtains  an  optimum  third-order  polynomial  least-squares  approximation 
to  an  arbitrary  zero-memory  nonlinear  transformation.  Thus,  linear 
systems  in  parallel  with  square-law  systems  and  cubic  systems  are  of 
major  importance  in  modeling  many  physical  problems. 

Sections  2  and  3  discuss  nonlinear  wave  force  problems  and  non¬ 
linear  drift  force  problems.  Appropriate  models  are  formulated  for 
each  problem  where  zero-memory  nonlinear  systems  are  followed  by 
linear  systems.  Two  types  of  problems  are  solved:  (a)  the  spectral 
decomposition  problem,  and  (b)  the  system  identification  problem. 
Statistical  errors  in  estimates  are  given  to  help  design  experiments 
and  to  evaluate  measured  results.  Section  4  develops  input/output 
relationships  for  stationary  random  data  through  three  zero-memory 
nonlinear  nonsymmetrical  systems:  (a)  three-slope  systems,  (b) 

catenary  systems  and  (c)  smooth-limiter  systems.  Formulas  are  stated 
for  input/output  probability  density  functions,  correlation  functions 
and  spectral  density  functions. 
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1.  PARALLEL  LINEAR  AND  NONLINEAR  SYSTEMS 


Before  discussing  techniques  for  analyzing  and  identifying  system 
properties  In  parallel  linear  and  nonlinear  systems,  three  topics  need 
to  be  covered  that  will  occur  In  later  material.  These  topics  are 

1.  zero-memory  nonlinear  systems, 

2.  optimum  third-order  polynomial  least-squares  approximations  to 
zero-memory  nonlinear  systems,  and 

3.  finite-memory  nonlinear  systems. 


1.1  ZERO-MEMORY  NONLINEAR  SYSTEMS 

A  "zero-memory"  nonlinear  system  Is  a  system  that  acts 
Instantaneously  on  present  Inputs  In  some  nonlinear  fashion.  It  does 
not  weight  past  Inputs  due  to  "memory"  operations  as  In  convolution 
integrals  for  constant  parameter  linear  systems,  where  the  present 
value  of  the  output  Is  a  function  of  both  present  and  past  values  of 
the  Input.  For  zero-memory  nonlinear  systems,  the  output  y(t)  at  any 
time  t  Is  a  single-valued  nonlinear  function  g[x(t)3  of  the  input  x(t) 
at  the  same  Instant  of  time.  Thus,  as  shown  In  Figure  1.1, 


yCO  s  *tx(t))  (i.i) 

where  g(x)  is  a  single-valued  nonlinear  function  of  x.  Note  that  g(x) 
is  not  a  function  of  t. 
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x(t) 


Nonlinear 

System 

g(x) 


y(t) 


g[x(t)] 


Figure  1.1  Zero  memory  nonlinear  system. 


The  system  g(x)  is  nonlinear  means  that  for  any  constants  a^,  a^ 
and  any  inputs  Xj,  Xj 

g(aJx1  ♦  a2x2)  t  *l  g(xl)  +  a2  g(x2)  (1.2) 

either  because  g(x)  is  not  additive  or  because  g(x)  is  not  homogenous,  or 
both.  The  system  is  a  constant  parameter  nonlinear  system  if  the  response 
of  the  system  is  independent  of  the  particular  time  of  use,  namely,  if 
an  input  x(t)  is  translated  to  an  input  x(t  ♦  t),  then  the  output  y(t)  is 
translated  to  an  output  y(t  ♦  t).  If  the  system  is  a  constant  parameter 
nonlinear  system,  and  if  the  input  x(t)  represents  a  stationary  random 
process,  then  the  output  VCO  will  also  be  a  stationary  random  process. 
For  these  cases  y(t)  *  glx(t)]  gives  y(t  ♦  x)  *  g[x(t  ♦  x))  and  the 
correlation  functions 

Rxy(T)  58  E[x(t)y(t  +  x )]  *  E{x(t)g[x(t)  +  x)]}  (1.3) 

Ryy(0  ■  E[y(t)y(t  +  t)]  *  E{g[x(t)]g[x(t  +  x )] }  (1.4) 

are  functions  only  of  t  as  required  for  stationary  random  data. 
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Many  examples  of  zero-memory  nonlinear  symmetrical  systems  are 
discussed  in  References  [1,  2].  They  include  two-slope  systems,  dead- 
zone  systems,  clipped  systems,  square-law  systems,  cubic  systems, 
hardening  spring  systems  and  softening  spring  systems.  Three  examples 
of  zero-memory  nonlinear  nonsymmetrical  systems  are  treated  here  in 
Section  4. 

Three  theorems  are  proved  in  References  [1,  2]  that  are  useful 
for  determining  input/output  relations  when  stationary  random  data 
pass  through  a  zero-memory  nonlinear  system.  Theorem  1  applies  to 
arbitrary  stationary  random  input  data.  Theorems  2  and  3  apply  only 
to  Gaussian  stationary  random  input  data.  These  theorems  will  be  used 
in  Section  4. 


Theorem  1. 


For  any  input  data  x(t)  with  probability  density  function  p(x)  where  x(t) 
passes  through  a  zero-memory  nonlinear  system  to  produce  y  =  g(x)  which 
is  single-valued  and  one-to-one,  the  output  probability  density  function 
P2(y)  for  the  output  y(t)  satisfies  the  relation 


p2(y) 


=  •  P(*> 
|dy/dx| 


(1.5) 


This  theorem  assumes  that  (dy/dx)  =  g'(x)  exists  and  is  not  equal  to 
zero.  When  solving  for  p2(y),  the  variable  x  on  the  right-hand  side 
should  be  replaced  by  its  equivalent  y  from  x  =  g~*(y). 
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Theoren  2.  (Price) 


For  Gaussian  input  data  x(t)  with  known  autocorrelation  function 
RXXCX).  where  x(t)  passes  through  any  zero-memory  nonlinear  system  to 
produce  y  =  g(x),  the  output  autocorrelation  function  R^Cx)  for  y(t) 
satisfies  the  relation 


3R, 

3R 


xx 


(X) 

TxT 


Efg'Cxpg'Cxj)] 


(1.6) 


whenever  g'(x)  -  [dg(x)/dx]  exists  at  xx  -  x(t)  and  x2  -  x(t  +  T)  . 

Theorem  3.  (Bussgang) 

For  Gaussian  input  data  x(t)  with  knhwn  autocorrelation  function 
ftxx(T)*  where  Passes  through  any  zero  memory  nonlinear  system  to 
produce  y  =  g(x),  the  input/output  cross-correlation  function  R  ft) 

XyV  W 

satisfies  the  relation 


Rxy(x)  #  ^  |  *  gCxJpCxJdx-^I^ErxgOt)]  (1.7) 

In  Theorems  2  and  3,  the  first-order  probability  density  function 
for  x  Is  given  by  the  Gaussian  form 

P(x)  ■  (l/o/27)exp(-x2/2o^  (1.8) 

where  the  mean  value  Is  assumed  to  be  zero  and  the  varl- 
2  7 

anCe  °x  =  E£xc(t)].  In  Equation  1.7,  the  quantity 

«l 

E[xg(x)]  *  /  xg(x)p{x)  dx  (1.9) 
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Identification  of  Zero-Memory  Nonlinear  Systems 


When  y  =  g(x)  is  single-valued  and  one-to-one,  it  can  be 
identified  from  measurements  of  x(t)  and  y(t)  by  using  the  relation 

y  x 

P2(y0)  s  /  p2(y)4y  =  /  p(x)dx  =  p(xq)  (i.io) 

-«0  —06 

where  P (x0 )  and  P2(y0)  are  the  Probat>i  1  i ty  distribution  functions  of 
x(t)  at  xQ  and  of  y(t)  at  yQ  =  g(xQ),  respectively.  To  determine  the 
zero-memory  nonlinear  function  y  =  g(x),  one  should  select  various 
values  of  xQ,  calculate  P (xQ)  and  then  determine  the  associated  values 
of  yQ  such  that  £(yQ)  =  P(xQ) . 


1.2  OPTIMUM  THIRD-ORDER  POLYNOMIAL  APPROXIMATION 

Suppose  y  =  g(x)  represents  any  zero-memory  nonlinear  system 
where  y  =  y(t)  and  x  =  x(t).  What  Is  the  optimum  least-squares 
approximation  to  y  =  g(x)  by  the  third-order  polynomial  y  =  y(t)  where 


y  * 


aix 


2  3 

+  a2x  +  a^x 


(1.11) 


under  the  assumption  that  x  follows  a  zero  mean  value  Gaussian 
distribution?  To  be  specific,  what  should  be  the  choices  of  the 
coefficients  a^,  a2  and  aj  so  as  to  minimize  the  quantity 

Q  =  E[y  -  y)2]  =  E[  (y  -  a^  -  a2x2  -  a3x3)2]  (1.12) 


over  all  possible  choices  of  these  coefficients? 
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For  any  y,  the  Gaussian  assumption  on  x  gives 


Q  =  E(y2)  -  2a.^E(xy]  -  2«^E[x2y)  -  2a3Elx3y] 


*  +  ♦  tt|E(x4J  +  a*E(x‘]  (1.13) 

using  the  fact  that  E[x]  *  E[x3]  -  E[x5]  -  0  for  the  p(x)  of  Equation 
1.8.  Then  setting  partial  derivatives  of  Q  with  respect  to  aj,  a2. 


and  a3  equal  to  zero  shows 

4^-  =  -2E[xy)  +  2axElx2)  ♦  2<r3Elx4)  a  o 

(l.lU) 

3  -2Elx2y)  ♦  ai^lx4)  a  0 

(1.15) 

|a*3  3  -22l*3yl  ♦  Sa^lx4)  ♦  Sa^l*6)  a  0 

(1.16) 

Now.  using  the  fact  that  E[x2]  -  o2,  Elx4]-  3oJ.  and  E[x6] 
the  Gaussian  p(x)  of  Equation  1.8  shows 

-  15ojj  for 

ai°i  +  Vi  *  EUyJ 

(1.17) 

3a2o4  a  E(x2y) 

(1.18) 

**A  *  15*fl  *  «U3yl 

(1.19) 
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Hence,  for  any  y,  one  obtains  the  y  of  Equation  1.11  from  the 
special  coefficients 


15o^EtxyJ  -  3o4E(x3y] 

(1.20) 

1  ~  « 

a2  = 

2  3®: 

(1.21) 

o2E[x3y]  -  3ojE[^y] 

a3  *  71 

^x 

(1.22) 

Example.  Square-Law  System  with  Sign 

Consider  application  of  these  matters  to  a  square-law 

systea  with 

sign  where 

y  ■  x|x| 

(1.23) 

For  this  particular  y,  the  Gaussian  assumption  on  x  yield 

ElxyJ  *  E[x2|x|]  *  2E[x3)x>0  *  2o3\|(2/n) 

(1.2**) 

E(x2y)  *  E(x3|x|]  *  o 

(1.25) 

E(x3yl  8  EIx4|x|J  *  2E[x5)x>q  s  8o3NJ(2/n) 

(1.26) 
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From  Equations  (1.20)  to  (1.22),  one  obtains 


al  “ 

(1.27) 

a2  *  0 

(1.28) 

a3  *  / ( 2/ir ) /3 crx 

(1.29) 

Hence,  the  optimum  third-order  polynomial 

to  x | x |  by  Equation  1.11  Is 

~  3 

y  =  a^x  +  a^x 

least-squares  approximation 

s  jox  \J(2/n)  j  x  +  ^J(2/n)^3oxj  x3 

(1.30) 

Thus  x [ x |  can  be  treated  as  a  combination  of  a  linear  system  in 
parallel  with  a  cubic  system.  Figure  1.2  shows  x | x |  compared  to 
the  y  of  Equation  1.30. 
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Nonlinear  Output  Value 


1.3  FINITE-MEMORY  NONLINEAR  SYSTEMS 


When  finite -memory  operations  are  desired  in  connection  with 
zero-memory  nonlinear  systems,  they  can  often  be  obtained  by  inserting 
a  constant-parameter  linear  system  before  and/or  after  the  zero-memory 
nonlinear  system.  Cases  where  the  linear  system  Is  after  the  zero- 
memory  nonlinear  system  represent  the  cases  of  greatest  interest  and 
are  the  only  cases  discussed  in  this  report.  Other  cases  where  the 
linear  system  precedes  the  zero-memory  nonlinear  system  are  discussed 
in  References  [1,  2],  Figure  1.3  shows  a  finite-memory  nonlinear 
system  with  a  linear  system  A(f)  that  is  after  the  zero-memory 
nonlinear  system  defined  by  g(x). 


x(t) 


Figure  1.3  Finite-memory  nonlinear  system. 
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An  extension  of  Figure  1.3  that  is  applicable  to  many  physical 
problems  occurs  when  the  input  data  x ( t )  passes  through  parallel 
linear  and  finite-memory  nonlinear  systems  as  shown  in  Figure  1.4, 
where  H(f)  and  A (f )  are  two  arbitrary  constant-parameter  linear 
systems  and  where  g(x)  is  any  specified  zero-memory  nonlinear 
system.  The  output  noise  data  n(t)  is  assumed  to  be  uncorrelated  with 
both  x(t)  and  y(t). 


x(t) 


y(t) 


Figure  1.4  Parallel  linear  and  finite-memory  nonlinear  systems. 
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Two  problems  are  associated  with  the  nonlinear  model  in  Figure 

1.4. 

(1)  Spectral  Decomposition  Problem.  Given  H(f),  A(f)  and  g(x),  plus 
measurements  only  of  x(t),  determine  the  spectral  properties  of 
y 1  (t )  and  y2(t).  If  y(t)  is  measured  as  well  as  x(t),  determine 
also  the  spectral  properties  of  n(t). 

(2)  System  Identification  Problem.  From  simultaneous  measurements  of 
both  x(t)  and  y(t),  identify  the  optimum  frequency  response 
functions  H(f)  and  A(f)  to  minimize  the  autospectrum  of  n(t). 

Sections  1.4  and  1.5  outline  practical  techniques  that  can  solve  these 
two  problems.  In  Section  1.4,  input  data  can  be  non-Gaussian,  whereas 
in  Section  1.5,  input  data  are  assumed  to  be  Gaussian.  Thus,  these 
techniques  are  applicable  to  physical  situations  of  widespread 
importance  without  restricting  the  probability  density  function  of  the 
measured  input  data.  More  general  models  than  Figure  1.4  are 
considered  involving  one  linear  path  and  two  parallel  nonlinear  paths 
as  shown  in  Figure  1.5.  Sections  2  and  3  apply  these  techniques  to 
wave  force  problems  and  drift  force  problems. 
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1.4  N0N-6AUSSI AN  INPUT  DATA 


Consider  the  general  single-input/single-output  nonlinear  model 
of  Figure  1.5  with  three  parallel  paths  where  the  input  data  x(t)  can 
be  non-Gaussian.  Let  Q2CX (t )3  be  an  arbitrary  known  zero-memory 
nonlinear  transformation  of  x (t )  and  let  gjCx (t ) 3  be  a  different 
arbitrary  known  zero-memory  nonlinear  transformation  of  x(t).  Let 

Xl(t)  =  x (t) ,  x  2 ( t )  =  g2[x(t)],  x  3 ( t )  =  g3[x(t)]  (1.31) 

represent  the  three  usually  correlated  input  records  to  the  three 
linear  systems  Aj(f),  A2(f)  and  A3 (f ) ,  respectively.  The  three 
associated  usually  correlated  output  records  from  these  systems  are 
denoted  by  y ^ (t ) .  y2(t)  and  y3(t ) ,  respectively.  To  complete  the 
model,  let  n(t)  represent  extraneous  uncorrelated  output  noise  and  let 
y(t)  represent  the  total  output  from  the  system.  A  special  important 
case  of  Figure  1.5  is  when  g2[x{t)]  *  x2(t)  and  g3[x(t)]  =  x 3 ( t ) . 
This  case  is  treated  in  Section  1.5. 
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X(t) 


y(t) 


Figure  1.5  General  single-input/single-output  nonlinear  model 
for  non-Gaussian  input  data  that  passes  through  a 
linear  system  in  parallel  with  two  finite-memory 
nonlinear  systems. 


Note  that  Figure  1.5  simplifies  to  Figure  1.4  when  g3[x(t)]  *  0 
so  that  all  results  obtained  here  apply  to  the  model  in  Figure  1.4  by 
merely  setting  X3(t)  and  all  subsequent  terms  computed  from  X3 ( t )  to 
zero. 


Figure  1.5  can  be  replaced  by  the  equivalent  three-input/single¬ 
output  linear  model  of  Figure  1.6  where  the  capital  letter  quantities 
are  Fourier  transforms  of  associated  small  letter  quantities.  To  be 
specific,  let 
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Xl(0  ■  Jfoltt)].  »i(f) 

x2(f)  ■  ^C*2(t)].  V2(f) 

»3(f)  *  5?C*3(‘)3.  Y3(f) 

H(f )  ■  5^[n(t)].  Y(f) 


(1.32) 

■  "pin  (t)] 

(1.33) 

•pi >3(t)] 

(1.34) 

•pty(t)} 

(1.35) 

Measurement  of  x(t)  and  y(t)  enables  one  to  compute  the  quantities 
xl(f)i  X2(f),  X3(f)  and  Y(Y)  when  x2(t)  and  X3 (t )  are  known. 

Recognition  of  the  equivalence  between  Figures  1.5  and  1.6  is  a 
significant  achievement  because  Figure  1.6  can  be  solved  by  well-known 
multiple-input/single-output  techniques  that  are  derived  and  discussed 
fully  in  References  [3,  4],  These  procedures  are  applicable  for  input 
data  that  can  be  Gaussian  or  non-Gausslan.  Independent  research  on 
these  matters  based  on  References  [5  -7]  is  in  Reference  [«. 


xx(f) 

X2(f) 

X3(f) 


Figure  1.6  Equivalent  three-input/single-output  linear  model 
to  Figure  1.5  where  the  three  input  records  can 
be  correlated. 
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The  basis  of  multiple-input/single-output  procedures  for  solving 
Figure  1.6  is  to  change  the  input  records  by  conditioned  spectral 
density  techniques  so  that  Xx(f)  is  left  alone,  X2(f)  is  changed  to 
X2.i(f)  where  the  linear  effects  of  X ^ (f }  are  removed  from  X2(f),  and 
X3 ( f )  is  changed  to  ^3 . 2 1  (f )  where  the  linear  effects  of  both  X^ (f ) 
and  X2(f)  are  removed  from  X3(f).  These  new  input  records,  Xj(f), 
X2.1(f)  and  x3- 2 ! (^ )  will  then  be  mutually  uncorrelated  and  become  the 
inputs  to  the  revised  three-input/single-output  linear  model  shown  in 
Figure  1.7.  The  noise  output  record  N(f)  and  the  total  signal  output 
record  Y ( f )  are  the  same  as  before.  However,  the  three  previous 
separate  output  records  Y ^ (f ) ,  Y2(f)  and  Y3 (f )  are  now  replaced  by 
three  new  separate  output  records  Ya(f),  Y^(f)  and  Yc(f)  that  will  be 

mutually  uncorrelated.  Also  the  three  previous  linear  systems  Aj(f), 

A2(f)  and  A3 (f )  are  now  replaced  by  three  new  linear  systems  Lj ( f ) , 
L2(f)  and  1.3(f). 


Figure  1.7  Revised  three-input/single-output  linear  model 
equivalent  to  Figure  1.5  where  the  input  records 
are  mutually  uncorrelated. 
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To  simplify  the  notation,  the  three  mutually  uncorrelated  input 
records  in  Figure  1.7  will  be  denoted  by  ll^(f),  ^(f)  and  113(f)  where 


Ui(f)  =  Xi(f)  «  X(f) 

(1.36) 

U2(f)  •  *2-l(f) 

(1.37) 

U3(f)  =  X3-2!(f) 

(1.38) 

Figure  1.7  is  now  essentially  three  separate  single-input/single- 
output  linear  models  where  the  linear  systems  Lj(f),  12(f)  and  1.3(f) 
can  be  computed  by  the  usual  spectral  relations 


Vf) 

Li(f)  *  dlTT 


uiui 


vf} 

Lz(n  *  — 


u2u2 


(f) 


G  (f) 
u3yv  1 

L*(f)  s  rSfT 

11  .  11  v  • 


U3U3 


(1.39) 


(1.40) 


(1.41) 


The  6(f)  quantities  are  one-sided  spectral  density  functions  that  can 
be  computed  easily  using  formulas  in  References  [3,  4]. 
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Note  that  the  linear  system  L^f)  is  the  same  as  the  overall  optimum 
linear  system  HQ(f)  between  the  input  x(t)  and  the  total  output  y(t) 
as  computed  by 


G  (f) 

Vf>’m  <‘-42> 

xxv  ' 

where  G ( f )  is  the  cross-spectral  density  function  between  x(t)  and 
y(t),  and  Gxx(f)  is  the  autospectral  density  function  of  x(t). 

Other  relations  for  Figure  1.7  are  as  follows. 


Y(f)  *  Ya(f)  +  Yb(f )  ♦  Yc(f )  ♦  N(f) 

(1.43) 

G«(f)  ■  V.(f)  +  V  +  SVc,f)  +  G""(f) 

(1.44) 

V.(f)  ■  l4(Ol\Ul(0  ”  r2iy(f)syy(f) 

(1.45) 

Vb(f)  =  lL2<f)|Zsu2u2<<')  ■  Yu2y(f)Gyy(f) 

(1.46) 

Vc(f)  *  lL3<f)l\u3(f)  * 

(1.47) 

Gnn<f>  ‘  I1  -  Y^tf)  -  Y^tf)  -  Y^O^f) 

(1.48) 

2  2  2 

The  quantities  y  ..(f),  (f)  and  y  (f)  are  the  ordinary  coherence 

V  u2y  u3y 

functions  between  y(t)  and  the  three  inputs  u^(t),  u2 (t )  and  U3U), 


respectively. 
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The  quantities  U2(f)  and  113(f)  in  Figure  1.7  are  given  by 


U2(f)  =  X2(f)  -  L^OX^f)  (1.49) 


U3(f)  =  X3(f)  -  L^mX^f)  -  L23(f)U2(f)  (1.50) 

where 


6l2(f) 

L12(f)  =  GJ7(T)  t1*51) 


613(f) 

Li3(f)=^TfT  (1-52) 


,  _  G23^)  *  Li3(f)G2iCf) 

L23(f)  "  622(f)  -  L12(f)G21(f) 


From  knowledge  of  L1(f),  L2(f)  and  L3(f),  the  linear  systems  Ax (f ) , 
A2(f)  and  A3(f)  can  be  computed  by  the  algebraic  relations 


A3(f)  ■  l3(f) 

(1.54) 

A2(f)  =  L2(f)  -  L23(f)A3(f) 

(1.55) 

AL(f)  =  L x ( f )  -  L12(f)A2(f)  -  L13(f)A3(f) 

(1.56) 
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In  terms  of  more  basic  spectral  quantities,  the  terms  used  in 


Equations  1.39  to  1.41  to  compute  L^(f)  L.2(f)  and  L3 (f )  are 


V<f)  ■  Gly(f) 

(1.57) 

Vi(f)  =  Gi‘(,) 

(1.58) 

Gu2y(f)  '  G2ylf)  '  Gl(f)G2l(f) 

(1.59) 

Gu2u2(f)  ’  G22(f>  -  >-i2(f)G21(f) 

(1.60) 

su3y(f )  -  G3y(f)  -  LjfflGjjlf)  -  L2(f)G3u2(f) 

(1.61) 

Gu3u3(f)  *  G33<f)  -  <-l3tf)G31{f)  -  l23(f>G3u2in 

(1.62) 

where 

G3u2(f)  *  632(f)  -  L12(f)G31(f) 

(1.63) 

All  of  these  terms  are  generally  required  when 

Gaussian  input  data.  Note  that  Gy  u  (f)  t  622(f)  and 

dealing  with  non- 

6u3u3(f>  *  G33(H 

except  for  special  cases. 
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Some  further  points  are  worthy  of  note.  The  relation 


4(f)  =  Ax(f)  +  L12(f)A2(f)  +  L13(f)A3(f)  (1.64) 

shows  that,  in  general,  Lj(f)  t  Aj(f).  Thus  determination  of  Li(f), 
representing  the  overall  optimum  linear  system  HQ ( f )  between  Xj(f)  and 
Y(f),  is  not  the  same  as  the  actual  linear  system  Aj(f)  that  exists 
between  X^(f)  and  Y^(f).  The  spectral  output  of  y^(t)  given  by 

Gy  y  (f)  =  |A1(f)|2Gu(f)  (1.65) 

will  generally  differ  from  the  spectral  output  of  ya(t)  given  by 

Gv  y  (,)  "  j*-i(f))2Gi1(f)  (1.66) 

Note  also  that  the  relation 

L2(f)  »  A2(f)  +  43(f)A3(f)  (1.67) 

shows  that  in  general  that  L2(f)  /  A2(f).  The  system  l2(f)  represents 
the  overall  optimum  linear  system  between  U2(f)  and  Y(f).  This  is  not 
the  same  as  the  actual  linear  system  A2(f)  that  exists  between  X2(f) 
and  Y2(f).  Also,  the  spectral  output  of  y2(t)  given  by 

sy  y  (f)  =  lA2(f>l2G22^f>  (1’68> 
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will  generally  differ  from  the  spectral  output  of  y^t)  given  by 


Gy  y  (f)  -  |L2(f)re  u  (f) 

ybyb  c  22 


(1.69) 


However,  note  the  equality  relation  that 


1.3(f)  -  A3(f) 


(1.70) 


where  13(f)  represents  the  overall  optimum  linear  system  between  U3(f) 
and  Y(f),  while  A3(f)  is  the  actual  linear  system  between  X3(f)  and 
Y3(f).  In  general,  since  G  (f)  /  G33(f),  the  spectral  output  of 

U3U3 

y3(t)  given  by 

(1.71) 


%y3<f>  =  lA3<f)l  G33<f) 


will  differ  from  the  spectral  output  of  yr(t)  given  by 


Vc'0  '  |L3(f)l  Su3<f) 


(1.72) 


Example  of  Figure  1.4 

Consider  the  example  of  Figure  1.4  where  a  linear  system  Is  in 
parallel  with  only  one  finite-memory  nonlinear  system.  This  example 
corresponds  to  Figure  1.5  where  H(f)  *  Aj(f),  g(x)  =  g2(x),  A(f)  = 
A2(f)  and  g3(x)  *  0.  Appropriate  formulas  for  this  example  are  as 
follows. 
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Ux(f)  =  Xj(f) 

(1.73) 

U2(f)  =  X2(f)  -  L12(f)X1(f) 

(1.74) 

l~l  2  C  ^ )  =  CGi2(f  )/6i^(f )] 

(1.75) 

4(f)  =  CGly(f)/Gn(f)] 

(1.76) 

• 

L2(f)  =  [SU2y(f)/GU2U2(f)J 

(1.77) 

where 

GU2y(f)  =  G2y(f)  -  L*2(f)Gly(f) 

(1.78) 

GU2u2(f)  -  G22(f)  -  Lu(f)G21(f) 

(1.79) 

Finally,  for  the  systems  in  Figure  1.3, 

A(f)  »  A2(f)  *  L2(f) 

(1.80) 

H(f)  =  Aj(f)  -  4(f)  -  l12(f)A(f) 

(1.81) 

Note  that  H(f )  /  4(f)  unless  A(f)  «  0.  Individual 

spectral  outputs 

are 

«,  <f>  *  |H(f)|ZG11(f) 

(1.82) 

Gy2y2(f)  *  |A(f)|2G22(f) 

(1.83) 

Gy  y  (f)  ■  l4(f)|2GU(f) 

(1.84) 

Vb(f)  '  |L*(f,lV2(f) 

(1.85) 

-25- 


The  total  output  spectral  density  function  is  given  by 


Gyy(f)  » 


Va 


(f)  + 


G  (f) 

Vb 


where 

•  GVa(f)  -  ,5iy(f)Gyy(f) 

Vb(f) '  v(f) vf) 

G„„<f>  *  U  -  V<fl  -  V<f>lGyy<f> 


(1.86) 

(1.87) 

(1.88) 

(1.89) 


These  formulas  are  general  results  for  Figure  1.4  where  the  input  data 
can  be  non-Gaussian  and  where  X£ (f )  represents  the  Fourier  transform 
of  any  zero-memory  nonlinear  system  output  x2(t)  *  g[x(t)]. 


This  concludes  the  example  of  Figure  1.4. 


1.5  GAUSSIAN  INPUT  DATA 


A  special  case  of  great  interest  is  the  case  where  input  data  are 
Gaussian  and  where  the  zero-memory  nonlinear  systems  in  Figure  1.4  are 

g2Cx(t)]  =  x2(t),  g3[x (t)]  =  x3(t)  (1.90) 

This  case,  shown  in  Figure  1.8,  is  the  Case  1  single-input/single¬ 
output  nonlinear  model  in  References  [l,  2]. 


y(t 


Figure  1.8  Case  1  single-input/single-output  nonlinear  model 
with  Gaussian  input  data. 


The  Fourier  transforms  X2(f)  and  X3 ( f )  on  the  related  three- 
input/single-output  linear  model  of  Figure  1.6  have  a  very  special 
meaning  for  this  Case  1  nonlinear  model.  Specifically,  X2 (f )  is  the 
Fourier  transform  of  x^(t)  and  X3 (f )  is  the  Fourier  transform  of  x3(t) 
as  denoted  by 

X2(f)  -  ?[x*(t)]  (1.91) 

X3(f)  ■  ?[x3(t)]  (1.92) 

For  Gaussian  input  data,  x(t)  will  be  uncorrelated  with  x^(t)  but  x(t) 
will  be  correlated  with  x3(t).  It  follows  that  the  cross-spectral 
density  functions  G^2(f)  =  0  and  G23(f)  =  A1  sc^  as  shown  in 

Reference  [1],  the  cross-spectral  density  function 

Gu(f)  =  3o^Gu(f)  (1.93) 

where  the  variance  o2  *  E[x2(t)]  when  E[x(t)|  *  0.  Here,  L12(f)  s  0 
and  L23(f)  *  0  but 

L13(f)  =  3o2  (1.94) 

The  three  mutually  uncorrelated  input  records  Uj(f),  U2(f)  and 
U3(f)  in  Figure  1.7  are  given  now  by  the  simple  relations 

Uj(f)  -  Xx(f)  (1.95) 

U2(f)  “  X2(f)  (1.96) 

U3(f)  «  X3(f)  -  ^(f)  (1.97) 

Thus,  Figure  1.7  becomes  Figure  1.9  for  the  model  of  Figure  1.8,  and 
the  L-systems  can  be  computed  easily  by  Equations  1.39  to  1.41. 
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Figure  1.9  Revised  three-lnput/slngle-output  linear  model 
equivalent  to  Figure  1-8  where  the  input  records 
are  mutually  uncorrelated. 


Specifically  for  the  model  of  Figure  1.9,  the  L-systems  are  given  by 
the  formulas 


where 


Ll(f) 

.  S(f) 

(1.98) 

.  Vf) 

(1.99) 

4(f) 

Vf) 

W* 

(1.100) 

V(f)  ■  V°  -  ^V0 

(1.101) 

U3u3^f^  *  G33^^  "  lCf ^ 

(1.102) 
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In  place  of  Equations  1.54  to  1.56,  simpler  relations  between  the 
linear  systems  in  Figures  1.8  and  1.9  are  now 


A3(f)  =  L3{f) 

(1.103) 

A2(f)  =  L2(f) 

(1.104) 

Aj ( f )  =  Ljff)  -  3o^A3(f) 

(1.105) 

Note  that  computation  of  Aj(f)  requires  knowledge  of  A3 ( f ), 
Al(f)  ^  4(f)  when  A3 ( f )  j4  0.  Also  the  relation 

and  that 

4(f)  *  Aj(f)  *  3»^A3(f) 

(1.106) 

2 

shows  that  L^(f)  is  a  function  of  the  input  variance  and  so  will 
change  with  different  Input  records.  However,  the  systems  Aj (f ), 
Ag ( f )  and  A3(f)  in  Figure  1.8  are  independent  of  the  input  variance 
when  Figure  1.8  Is  a  valid  nonlinear  model.  Thus,  conventional  linear 
system  identification  techniques  where  only  Lj(f)  is  computed  give 
erroneous  estimates  of  Aj(f). 

This  concludes  Section  1. 
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2.  NONLINEAR  WAVE  FORCE  MODELS 


Section  2  details  a  methodology  for  analyzing  parallel  linear 
and  nonlinear  systems  involving  a  square-law  operation  with 
sign.  The  analysis  is  applied  specifically  to  the  problem  of 
decomposing  random  wave  forces  on  a  structure  into  linear 
(inertial  force)  and  nonlinear  (drag  force)  components.  A 
procedure  is  presented  for  identifying  the  individual  initial 
and  drag  force  parameters  based  solely  upon  measurements  of  the 
input  wave  velocity  and  the  output  wave  force.  The  input  wave 
velocity  is  assumed  to  be  a  Gaussian  stationary  random  process 
with  arbitrary  autospectral  density  function.  Formulas  are 
stated  for  random  errors  in  estimates.  These  formulas  are  also 
useful  in  the  design  of  experiments. 
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2.1  FORMULATION  OF  WAVE  FORCE  MODEL 


Morison's  Equation  is  widely  used  to  study  wave  forces  on 
structures,  References  [7,  where  a  wave  force  output  y(t)  due 

to  a  wave  velocity  input  x(t)  consists  of  two  main  components: 

(1)  an  inertial  force  m(t)  that  is  proportional  to  the  derivative 
x(t)  of  the  wave  velocity  x(t), 

(2)  a  ^rag_force  d(t)  that  is  proportional  to  x(t)|x(t)|. 

Morison's  Equation  is  thus  of  the  form 

y(t)  =  m(t)  +  d(t)  =  Cj  i(t)  +  C2  x ( t ) | x ( t ) |  (2.1) 

where  the  exact  mature  of  the  constants  Cj  and  C2  would  either  be 
known  or  to  be  determined.  In  Morison's  work,  Cj  and  C2  are  not 
functions  of  frequency.  It  is  clear  that  Figure  1.4  can  represent 
Equation  2.1  by  setting 


H(f)  -  C^TTf) 
A(f)  -  C2 


(2.2) 


Equation  2.1  can  apply  to  more  general  situations  by  allowing  and  to 
be  functions  of  frequency.  Figure  2.1  illustrates  this  wave  force  problem 
for  a  tall  structure  in  the  ocean. 
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The  problem  of  concern  is  to  determine  the  spectral  properties  of  the 
inertial  component  m(t)  and  the  drag  component  d(t)  from  measurements  of 
the  wave  velocity  x(t)  and  the  net  force  y(t) .  This  determination  is 
difficult  for  two  reasons: 


(a)  the  drag  force  d(t)  has  a  nonlinear  relationship  of  the  type 
x(t)|x(t)|,  namely,  a  square-law  operation  with  sign,  and 

(b)  the  inertial  force  m(t)  and  the  drag  force  d(t)  are  correlated. 


The  specific  model  to  be  evaluated  here  is  shown  in  Figure  2.2  where 

n(t)  is  extraneous  noise,  and  H(f)  and  A(f)  are  the  frequency  response 

functions  of  constant  parameter  linear  systems  given  by  the  Fourier  trans¬ 
forms  of  weighting  functions  h(r)  and  afr),  namely, 


H(f)  »7[h(T)]  -  f  h(Oe'j2nfTdx 

j°  (2.3) 

A(f)  ** 3^[aCT) ]  »  f  a(t)e  ^2flfxdT 


The  system  xjx|  is  a  zero  memory  nonlinear  system  with  an  instantaneous 
output  w(t)  *  x(t)  |x(t)  | .  Note  that  unlike  Morison's  relationship  in 
Equation  2.1,  the  systems  HCf)  and  A(f)  are  not  assumed  to  be  constants 
or  even  real  numbers. 
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Two  problems  will  be  treated  for  the  model  of  Figure  2.2. 

(1)  Spectral  Decomposition  Problem.  Given  H(f)  and  A(f)  plus. 

measurement  of  x(t)  only,  determine  the  spectral  properties  of  m(t) 
and  d(t).  If  y(t)  is  measured  as  well  as  x(t),  determine  also  the 
spectral  properties  of  n(t). 

(2)  System  Identification  Problem.  From  measurements  of  x(t)  and  y(t), 
identify  the  optimum  frequency  properties  of  H(f)  and  A(f)  to 
minimize  the  autospectrum  of  n(t). 

2.2  SPECTRAL  DECOMPOSITION  PROBLEM  ’ 

Referring  to  Figure  2.2,  assume  the  input  record  x(t)  is  a 
sample  time  history  of  a  stationary  (ergodic)  Gaussian  random 
process  with  mean  value  yx  =»  E[x(t)]  =  0.  Nonzero  mean  value 
inputs  are  treated  later  in  Section  2.U.  To  solve  the  Decomposition 
Problem  for  the  model  in  Figure  2.2,  an  approximation  is  required 
for  the  nonlinear  operation  x  |xj  .  The  third-order  polynomial  least- 
squares  approximation  is  given  by  using  Equation  1.30,  namely, 

xixi  ■  (°x.'il)x  +  (55;  '/D *3 

*  k(3oxx  +  x3)  k  ■ 

where  k  is  the  same  as  in  Equation  1.29  and  ox  is  the  standard 
deviation  of  x(t)  defined  by 


In  words,  x| x |  can  be  replaced  by  the  sum  of  a  linear  operation  plus  a 
cubic  operation.  For  an  input  x(t)  with  a  zero  mean  value,  this 
nonlinear  system  approximation  i s  Figure  2.3.  Substitution  of. 
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xft) 


Figure  2.3  Third-order  polynomial  for  x(t)  |  x(t)  |  . 


Figure  2.1*  Nonlinear  wave  force  model  with  correlated 
outputs  representing  Figure  2.2. 
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Figure  2.3  for  the  x|x|  operator  into  Figure  2.2  yields  the  model  shown 
in  Figure  2.**  which,  henceforth,  is  used  to  approximate  the  model  in  Figure  2*2« 
Formulas  will  now  be  derived  to  compute  the  autospectral  density  functions  of 
mCt)  and  d(t)  from  knowledge  of  H(f)  and  A(f) ,  using  measurements  of 
x(t)  and  y(t) . 

It  should  be  noted  here  that  Figure  2.1+  is  a  special  case  of  Figures  1.5 

and  1.8  that  can  be  analyzed  by  formulas  in  Sections  1.4  and  1.5.  However,  it 
is  more  instructive  to  directly  develop  here  the  desired  formulas  for  this 
special  problem.  In  agreement  with  physically  measurable  results,  one-sided 
spectral  density  functions  will  be  used  instead  of  theoretical  two-sided 
spectral  density  functions. 

In  Figure  2.4,  for  a  Gaussian  input  x(t)  with  zero  mean  value  and 
standard  deviation  cx,  the  terms  m(t) ,  d(t),  n(t)  and  y(t)  will  also 
have  zero  mean  values.  Formulas  for  the  Fourier  transforms  of  such  records 
with  long,  but  finite  length  T  are  given  by 


X(f)  =J*[x(t)]  =  (  x(t)e~^27rftdt 

;°  (2.6) 

V(f)  [y(tj]  -  [T  y(t)e“j2,rftdt 

Jo 


It  is  assumed  that  x(t)  and  y(t)  can  be  divided  into  n^  associated  sub¬ 
records,  each  of  length  T,  to  compute  desired  averages.  Autospectral  and 
cross -spectral  density  functions  of  x(t)  and  y(t)  are  then  computed  for 
one-sided  spectra  by 

Gxx(f)  -  |  e[x*  (f)X(f)] 

Gyy(«>  -  I  e[V  (f)Y(f)] 

Gxy(f)  -  |  E[x*(f)Y(f)] 


(2.7) 


Where  the  asterisk  (*)  denotes  complex  conjugate  end  E[  ]  is 

approximated  in  practice  by  an  ensemble  average  over  n 
quantities .  d 

m  terms  of  the  functions  H(f>  and  Mf>.  Fourier  transforms 

of  the  output  quantities  over  long,  but  finite  length  records 
yield 


V(f)  *  K<£)  +  D  { f )  +  N 
MCf)  =  H  { f )  x  ( f ) 


A(f)W(f) 


(2.8) 


(2.9) 


where  W(f)  is  the  Fourier  transform  o£  w(t)  de£.ned  by 

W(*>  =  k  X(£)  +  X3(fj]  (2 

Here,  w(t)  .  kC3ox2x(t)  +  x3(t)j  ^ 

x3«f>  =?[=<3<t)3  =  rx3(t)e-J2,,£tdt  (2.n) 

The  quantity  N(f)  is  the  Fourier  transform  „  .. 

r  transform  of  the-  unmeasured 

extraneous  output  noise  n(t)  that  is  assumed  to  be  uncorrelated  with 
both  m(t)  and  d(t). 

Figure  2.1,  Is  equivalent  to  the  tvo-input/single-output  linear  model 

shown  in  Figure  2.5  where  the  Inputs  x(t)  and  v(t)  are  correlated.  Figure 
2.5  is  a  special  case  of  Figure  1.6. 


(2.10) 
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Figure  2.5 


IlS"lnp!i/!t"9le'output  linear  mode1  with  correlated 
inputs  that  is  equivalent  to  Figure  2.4. 


cross-correlation 


Let  y(t>.  x3(t)  =  x3(t)  with  E[x(t) ]  «  0.  The 
function 


^xy(T)  =  e[x  (t)  y  {t+T)J  =  e[x  (t)x3  (t+T)j 

*  E[x(t)x(t+T)x(t+T)x{t  +  T)] 


For  zero  mean  value  Gaussian  data,  E[y(t)J  -  0  and  a  fourth-order 

moment  breaks  up  into  the  product  of  all  possible  pairs  of  second-order 
moments.  Specifically, 

eKw]  ■  «CvJ  eI>3-4]  *  «CvJ.  •[•*]♦ 

Hence,  R  (x  )  becomes 
xy 


Rxy(T)  -  3E[x(t+T)x(t+T)]  E[x(t)x(t+T)]  .  30^R  (T ) 
Th6  corresponding  cross-spectrum  relationship  is 


°xy(f)  -  3oxGxx<fl 


This  proves  by  letting.  y(t)  -  x  3(t)  »  x  3(t) 


that 


G  (f)  -  3crJo  ff) 


X  XX 


(2.12) 


Note  that 


W*'  ~  Gxx3(f)  -  Gx3x<£) 
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The  one-sided  autospectral  density  function  of  x^(t)  *  x  (t)  is 


Gx3X3Cf)  -  I  e[x3*(£)X3(£51 


(2.13) 


This  computation  can  be  done  directly.  The  one-sided  autospectral  dei 
function  of  w(t)  is  defined  by 


°ww(f)  -  I  e[w*(£)W<«>] 


(2.14) 


This  computation  requires  some  extra  work. 

A  formula  will  now  be  derived  for  Gww(f)  in  terms  of 

Gxx(f)  and  Gx  x  (f).  Equation  2.12  proves  that  x(t) 

3  3 

and  x^(t)  are  correlated  with 

gxx3<£>" 


(2.15) 


Hence, 


Gww <£>  -  kJl90iGxx(f>  +  30xgxx3(£>  ♦  3»xGxx3(£>  +  Gx3x3'£> 


‘H 

[21< 


>] 


G  _  ( f  )  +  G  (f ) 

X  XX  X jX J 


i 


(2.16) 


It  follows  also  that 


[• 


Gxw‘£>  -  M3<Gxx<£>  +  Gxx 


.">] 


-  6kc  G  (f) 
x  xx 


,(2.17) 


This  shows  that  x(t)  and  w(t)  are  correlated. 
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It  is  now  simple  to  compute  the  desired  autospectral  density 
functions  of  m(t)  and  d(t).  From  Equation  2.9» 

2 

Gmm  (f>  "  ^H(f)  I  Gxx(f)  (2’18) 

This  represents  the  inertial  term  in  Equation  2.1,  Also 

2 

Gdd(f)  -  |A(f)|  Gw(f)  (2.19) 

where  Gw(f)  should  be  computed  using  Equation  2.l6.  This 
represents  the  drag  term  in  Equation  2.1.  The  cross-spectral 
density  function  between  m (t)  and  d(t)  is  obtained  using 
Equation  2.17  as  follows. 

Gmd(f)  "  I  E[Ym*(f)Yd  {£)]  “  H*(f)A{f)  Gxw(f)  {2-20) 

»  6k H*  (f)Alf)Gxx.(f) 

Thus,  m  (t)  and  d  ( t)  are  correlated  and  from  Equation  2.8 


Cyy(f) 


nro 


(f)  ♦  Gdd(f)  ♦  GL,(f)  ♦  G^Cf)  ♦  Gnn(f) 


md 


mdv 


nn' 


, (2.21) 


Note 'that  computation  of  Equations  2.18  and  2.19,  involves  measurement  . 
only  of  x(t),  whereas  Equation  2.21  requires  measurement  also  of  y(t).  When 
both  x(t)  and  y(t)  are  measured,  Equation  2.21  can  compute  the  autospectrum 
of  the  unmeasured  n(t).  Good  models  occur  when  G  (f)  is  small  compared  to 
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From  Equations  2.1 6  and  2.19  the  drag  term  can  be  written 
Gdd(£)  -  k2|»(l)|2[^<Axlfl  +  “kjXj'*’]  (2 

This  (f )  is  the  sum  of  two  parts  obtained  by  using  the 

correlated  terms  x{t)  and  x^(t) .  An  alternative  formula  for 

G  (f)  is  derived  later  in  Equation  2.1*1*  where  G,.  (f) 
dd  uu 

is  based  upon  two  terms  that  are  uncorrelated.  From  Equation 
2.. 21,  note  that 

Gyy(f>  *  G™‘«  *  Gdd(»  *  Gnn(£)  ' 


(2.2U) 


2.3  SYSTEM  IDENTIFICATION  PROBLEM 


To  obtain  Gnn(f),  Equation  2.21  shows  that 

Gnn(£>  "Gyy(£)  -[Gm(fl  +  Gdd tf)  +  G*d  *  GId  ] 

A  more  informative  formula  for  Gnn(f)  will  be  derived  later 
in  Equation  2.1*9. 


The  starting  point  for  this  problem  is  the  model  in  Figure  2.1* 
where  it  is  now  assumed  the  properties  of  H(f)  and  A(f)  are  not 
known.  These  properties  will  be  determined  from  x(t)  and  y(t) 
based  upon  minimizing  the  spectral  density  function  Gnn(f)  of 
n(t)  over  all  possible  choices  of  linear  systems  to  predict  y(t) 
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from  x(t).  It  is  assumed  here  that  the  mean  value 

Px  **  Etx(t)D  “  0.  Nonzero  mean  value  inputs  are  treated  in 
Section  2.3. 

In  Figure  2.U,it  follows  that 


Y(f)  «  H  (f )  X  (f )  +  kA(f)  3o*X(f)  +  X3(f)  +  N(f) 


[3ohx 


(2.25) 


and  using  Equation  2.15, 

Gxy(f)  =  H(f)Gxx(f)  +  3k0xA(f  )Gxx(f)  +  kA  (f )  G  (f)  (2.26) 

xx3 

=  [h(£)  +  6ko3A(£)]  Gxx<£) 

Hence,  a,  proved  in  Reference  C23,  the  optimum  linear  system  is 
given  by 

G  ( f ) 

Vf)  “  5(77“  H(f)  +  6Ic^A(f)  (2-2T) 

xx  x 

This  system  HQ(f)  computed  from  x(t)  and  y(t)  only,  gives  the 
minimum' Gnn(f)  over  all  possible  linear  systems.  it  also  auto¬ 
matically  makes  n(t)  uncorrelated  with  x(t). 

Note  that  H  (f)  ^  H(f)  and  that  H  (f)  Is  a  function  of  the  input  variance 
o  o 

2 

(j  .  Thus  H  (f)  will  change  with  different  inputs  while  H(f)  will  be  the  same, 
x  o 

Note  also  that  determination  of  H(f)  requires  knowledge  of  A(f)  as  well  as 

H  (f).  By  using  H  (f)  instead  of  H(f),  Figure  2.U  can  be  redrawn  as  Figure 
o  o 

2.6.  This  is  then  equivalent  to  the  two-Input/single-output  linear  model 
of  Figure  2.7  where  x(t)  and  u(t)  are  uncorrelated.  Figure  2.7  is  a  special 
case  of  Figure  1.9.  Figure  2.7  can  also  be  derived  directly  from  Figure  2.5. 
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Fourier  transform  relations  for  Figures  2.6  and  2.7  show  that  Y(f) 

is  the  sum  of  a  linear  term  Y  (f),  an  uncorrelated  nonlinear  term  Y  (f), 

o  u 

and  an  uncorrelated  noise  term  N(f)  where 

yf>  "  H0(f)X(f)  -  j^H  { f )  +  6kojA(f)]  X(f)  (2.28) 

Yu(f)  "  A(f)U(f)  -  JcA(f)  [x3(f)  -  3ajjx(f)]  (2.29) 

Here,  the  cross-spectral  density  function 

G  (f)  «  0  (2.30) 

yo  yu 

Now,  noting  that  the  total  force  is  given  by 


Y(f)  «  Yq  (f)  +  Yu(f)  +  N(f)  (2.31) 

it  follows  from  Equation  2.30  that  the  autospectrum  of  the  total 
force  decomposes  into  three  additive  terms,  as  follows. 


GYY(f) 


y0y0 


(f)  +  G 


(f) 


+  Gnn(f> 


12.32) 


This  gives  the  model  in  Figure  2.6  that  is  equivalent  to  Figure  2.1*. 

The  quantity  u(t)  in  Figure  2.6  is  not  the  same  as  tfie  quantity 
w(t)  in  Figure  2.1*  and  U(f)>  W(f).  In  thia  cas6/ 

u(t)  »  k[x3(t)  -  3o*x(t)]  (2.33) 

which  has  a  Fourier  transform  from  a  long  but  finite  record  given 
by 


U(f) 


=  k 


(f) 


3o*X(f) 


(2.31*) 
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Here,  x(t)  and  u(t)  are  uncorrelated.  Using  Equation  2.15,  the 
autospectrum  of  u(t)  is  given  by 


^u<f)  -  **[oX3„ 

L  X3X 


(f)  "  3<Gvv.(f)  -  3o!g*„  (f) 

‘3 


3^3  x_xx3  -Vxx,'A'  +  9°xGxx(f) 


] 


3(f)  “  9°xGxx(f> 


(2.35) 


Also,  from  Equation  2.3U, 


Guy(f)  =  k[Gx3y(f)  "  3axGxy(f)]  (2-36) 

However, 

Guy(f)  =  A(f)Guu(f)  (2.3T) 


Hence,  it  follows  that 


A(f) 


Guy  (f ) 
Su  (f )  “ 


%Y<fl  -  3«xGxy») 

k[Gx3x3(f>  -  9»xGxx<£»] 


(2.38) 


where  computations  for  G  (f)  and  G  (f)  can  be  done  directly  from 

X3  3  x 3 y 

X  (f)  and  Y(f),  similar  to  direct  computations  of  G  (f)  and  G  (f)  from 
J  xx  xy 

X(f)  and  Y(f).  Equation  2.38  is  the  desired  result  to  identify  the  system 
A(f). 
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(2.39) 


The  system  H(f)  can  now  be  determined.  From  Equation  2.27,  using 
the  just  computed  A(f), 

H  (f )  =  H0(f)  -  6ko*A(f) 

This  is  the  desired  result  to  identify  the  system  H(f). 

After  A(f)  and  H(f)  are  computed,'  then  all  formulas  in 
Section  2.2  follow.  In  particular,  one  obtains  desired  inertial 
and  drag  results  by  using  Equations  2.1b  and  2.19.  From 
Equations  2.16  and  2.35, 

Gwv,<£>  -  36k2°xGxx<£)  +  Guu<£>  (£ 

Hence,  in  place  of  Equation  2.19,  one  can  write 


Su(£> 


G<uVf> 


+  CJ  .  .  .  (f) 
d  n 


(2.1*1) 


where  G  ^  ^  (f)  represents  the  linear  part  of  the  drag  term 

given  by 


"  36k^|A(f)|2Gxx(£) 


(2.1*2) 


and  G  ,  (f)  represents  the  uncorrelated  nonlinear  part  of 

n  d  n 

the  drag  term  given  by 


G.  ,  (f> 

d  n  ^  n 


|A(f)  TGUU  (f) 


(2.1*3) 


Here,  Guu(f)  should  be  computed  using  Equation  2.35.  Thus, 
G^ff)  is  the  sum  of  two  parts  obtained  from  the 
uncorrelated  terms  x(t)  and  u(t),  namely. 


<V£) 


A(f)  | 


36kMGxx|£> 


+  G 


uu 


(f) 


(2.1*1*) 
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E“  ^  iS  ^-tic  *.*.  further  neu 

results  for  the  eutospectre!  density  functions  of  yQ (t>  end 

(t)  in  Figures  2.6  and  2.7  are 

Gyoio(f)  *  !h0  <f)  |2Gxx(f}  (2.2*5) 

using  H0(f)  from  equation  2.67,  and 

Gyuyu(f)  =  iA(f,J\u(f)  (2.2*6) 

^sing  A(f)  from  Equation  2.38  and  G  (f|  f 

uu(f)  from  Equation  2.35. 

Note  that  Equation  o  2*6  { e  4.y,_ 

•  o  is  the  same  as  Equation  2.2*3. 

THe  linear  conference  function  y l  ,fl  „  „ 

xy(f)  between  the  input 
x(t)  and  the  total  output  y(t)  is  defined  by 


2  Ve‘f>  |Gxy<£>!2 

YXy  (f  )  =  — -  a  . . 

Gyy(£)  Gxx(£>Gyy<*> 


(2.2*7) 


where  the  G„  < 

y0yo(f)  is  9i^en  by  Equation  2.1*5.  nonl. 

_  .  o  ir,e  “Onimear 

coherence  function  q  (r)  <  c  ^  - . 

Y1'  is  defined  here  by 


G  (f ) 

y  y  V*-; 

a  •'u _ 

G  (fT 

yy v  ’ 


(2.2*8) 


where  the  Gy^yJf)  is  given  by  Eguation  2>Ji6>  _  Equation  2M  -s  the 

same  as  the  linear  coherence  function  between  the  uncorrelated  input  u(t)  and 
the  total  output  y(t). 


In  terms  of  these  coherence  functions,  the  output  noise 
autospectral  density  function  is  given  by  the  simple  formula 


G  (f) 
nn 


[>-' 


<y<£>  - 


2 

a  (f) 
xy v 


-] 


Gyy<f) 


(2.1*9) 


Clearly,  this  model  will  be  valid  at  frequencies  where  the  sum 
2  2 

tYxy(f)  +  ^^f)]  is  close  to  unity. 

2.1*  NONZERO  MEAN  VALUE  INPUT 

Suppose  that  the  input  x(t)  has  a  nonzero  mean  value 
UX  =  E[x (t )  ]  0.  For  the  case  where  x(t)  represents  the 

velocity  of  an  ocean  wave,  yx  would  be  the  underlying  current. 

When.  yx  £  0,  the- input  x(t)  should  be  expressed  as 

x(t)  -  [x(fc>  -  »x]  ♦  ux  <2-50) 

consisting  of  a  variable  input  term  [x(t)  -  Ux]  with  zero  mean 
value,  plus  a  constant  term  yx.  Now  the  cubing  operation 
x3(t)  becomes 

x3(t)  -  {[x(t>  -  yx]  +  yx)3  (2.51) 

-  [x(t)  -  uj  +  3(ux)[x(t)  -  „x]  +  3(ux)J[x(t)  -  ux].  (ux)3 

3 

where  the  final  constant  term  y  produces  an  effect  in  the 
spectrum  of  x3(t)  only  at  f  *  0«  For  f  £  0,  Equation  2. 51- 
shows  that  the  diagram  in  Figure  2.8  should  represent  the  nonlinear 
operation  x  (t)  for  the  input  [x(t)  -  Ux]«  Note  that  Figure 
2.8 simplifies  to  the  cubic  operator  by  itself  when  vx  *  0. 
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[x(t)-Mx] 


Figure  2.8  Nonlinear  system  for  x^(t)  when  ju  ^  0. 


Figure  2.9  Nonlinear  system  approximatlng.x(t) |x(t) 


when.  U  t  0 
x 
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Consider  next  the  nonlinear  operation  x(t)Jx(t)|  of  a 
square-law  system  with  sign  which  is  represented  by  the  • 
least-squares  approximation  of  Equation  2.1  ,  namely 

x(t)  |  x  ( t )  |  -  3ko^x  (t)  +  kx3(t) 


k  = 


3  0. 


4 


(2.52) 


(2.53) 


Here  from  Equation  2.50, 

x(t)  |  x(t)|  =  3koJ  {[x(t)-Mx]  +  Ux}  +  k{[x(t)-ux]  +  Mx 
-  k  [x(t)-Ux]  +  3k„x[x(t)-px] 

+  3k  (ojj  +  uj|>[x(t)-gx]  +  k(3o*ux+y^ 

where  the  final  constant  term  k(3o^yx  +  y^3)  produces  an 
effect  in  the  spectrum  only  at  f  *  0.  For  f  ^  0,  Equation  2.52 
shows  that  the  diagram  of  Figure  2.9  should  represent  the  nonlinear 
operation  x(t)  |x(t)j  for  the  input  Cx(t)  -  Ux3*  Note  that 
Figure  2.9  simplifies  to  Figure  2.3  when  =  0. 
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Return  now  to  the  original  problem  covered  by  Figure  2,k.  For 
situations  where  the  input  changes  to  [x(t)-y^]  with  yx  ^  0, 
Figure  2.k  becomes  the  new  extended  Figure  2.10.  For  Gaussian  data,  a 
linear  output  is  uncorrelated  with  the  output  from  a  square -law  operation 
but  is  correlated  with  the  output  from  a  cubic-law  operation.  Hence,  in 
place  of  Figure  2.6,  when  the  input  changes  to  [x(t)  -  yx]  with  ux  t  0, 
Figure  2.6  becomes  the  new  extended  Figure  2.11. 


The  optimum  linear  system  HQ(f)  in  Figure  2. 11  is  computed  as 
before  by  the  equation 


GXY(f) 

H  (f)  =  ~y 

°  Gxx(f) 

(2.5V) 

In  terms  of 

H ( f )  and  A(f),  one  finds  that 

H  (f)  =  H  ( f )  +  6ka*A(f)  +  3kU^A(f) 
o 

(2.55) 

•Note  this  reduces  to  Equation  2.27  when  Vx  *  0.  In 

Figure  2.11, 

the  Fourier 

transform  U(f)  of  u(t)  is  given  in  place 

of  Equation 

2.3^  by 

• 

V(f)  -  k[x3(f)  +  3yxX2(f)  -  3cxX(f)J 

(2.56) 

where 

X(f)  -7j[x(t)-Ux]} 

(2.57) 

X2(f>  ■=?■([  x(t)-ux]  } 

(2.58) 

x3(f)  ■^([x(t)-’Jx]  } 

(2.59) 
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One  should  now  proceed  as  wee  done  previously  to  compute 

G  (f)  and  G  (f)  where 
uu  uy 

Guu  (f)  =  |E[  U(f)  U(f)J  (2.60) 

Guy  (f)  -  |e[  U  { f )  y  ( f ) J  (2.61) 

Then  »  the  system  A(f)  is  computed  by 


A(f) 


Guy  <f> 
Guu  (f) 


and  from  Equation  2.55»  the  system  H(f)  is  computed  by 


(2.62) 


H(f)  88  H0(f)  -  6ko2A(f)  -  3kp2A{f)  (2.63) 

"  x 

The  inertial  component  m(t)  in  Figure  2.10  has  the  autospectrum 

Gm  {f)  85  I H  (f  >  J  2Gxx  (f )  (2.6U) 

The  drag  component  d  (t)  in  Figure  2.10  has  the  autospectrum 

Gdd<f>  “  |A(f)!2Gww(f)  (2.65) 


where  G  (f)  in  Figure  2.10  should  not  be  confused  with  G  (f)  in  Figure 
ww  uu 

2.11.  The  output  v(t)  in  Figure  2.10  is  quite  different  from  the  output 
u(t)  in  Figure  2.11.  Also,  the  output  m(t)  in  Figure  2.10  is  quite 
different  from  the  output  y  (t)  in  Figure  2.11. 
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The  quantity  G  (f)  is  computed  by  the  usual  formula 
w 


G  (f) 
w w 


ijE  [w*  ( f )  W  ( f  )J 


(2.66) 


where  W(f )  is  given  by 


W(f)  »  k  £x3  (f )  +  3yxX2<f)  +  3  (c?^+u^)X(f  )J 


(2.67) 


Note  that  Equation  2.67  reduces  to  Equation  2.10  when  px  a  o. 
The  relationship  between  U(f)  and  W(f)  is  given  by 


W(f)  -  U(f)  +  3k  (2o^+  M^)  X  (f ) 
which  leads  to  the  autospectral  density  function 


(2.68) 


6W(*)  *  =uu<f)  +  Gxx(f> 


(2.69) 


This  result  extends  Equation  2.40  for  yx  ^  0  and  occurs 


because 


G  (f)  =  0 
X2 


(2.70) 


Gxx3(£>  -  KGxx(f> 


Thus,  in  place  of  Equation  2.65,  one  can  write 


(2.71) 


G..(f)  »  |A(f)  |2G,iU  (f)  +  9k2(2o*+vj|)  lA(f)|4G  (f)  (2.72) 


X  X 


These  two  parts  of  G_.(f)  represent  uncorrelated  nonlinear 

QQ 

and  linear  parts  of  the  drag  term.  Equation  2.72  reduces  to 
Equation  2.44  when  yx  *  0. 
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Formulas  derived  in  Section  2.4forPx  /o  have  been  shown 
to  be  extensions  of  corresponding  formulas  in  Section  2.2  and  2.3 

vhere  Px  »  0.  In  place  of  Equations  2. 1*5  and  2.1*6,  when 

ux  t  0, 


Gy  y  (f>  -  lHo  <f > 1  Gv*<f> 


(2.73) 


using  HQ(f)  from  Equation  2.51*,  and 


o  (f)  =  |a(£)|2Guu(£) 

•'u  JU 


(2.7U) 


using  A(£)  £rom  Equation  2.62  and  0j]u(f)  trom.Equation  2  6o_ 
The  linear  coherence  function  Y^y(£)  ia  here 


,  Gv  v(f) 

Y2  If)  =  °\ 

*y(£)  Gyy(fj 


(2.75) 


using  the  Gy^ff)  from  Equation  2.73,  and  the  nonlinear 

■  e  -  O  .  . 


coherence  function  q2  (f)  is 

xy 

v(f) 


Gy  y  (f) 

Gyy(f> 


(2.76) 


using  the  Gy^yjf)  from  Equation  2.7k.  m  terms  of  the 
coherence  functions  from  Equations  2.75  and  2.76,  the  output 
noise  autospectral  density  function  is  now 


Gnn(£)  ’ 


(2.77) 
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2.5  RANDOM  ERRORS  IN  ESTIMATES 


Referring  back  to  Figure  2.5  and  Equations  2.68  through  2.32, 
one  can  view  the  model  in  Figure  2.5  as  a  two-input  system  with 
uncorrelated  inputs  as  shown  in  Figure  2.7. In  this  form,  the 
random  errors  to  be  expected  in  the  computation  'of  certain 
critical  parameter  values  are  given  directly  from  Reference [1] 
In  terms  of  the  normalized  random  error  (coefficient  of 

A 

variation)  of  an  estimate  <J>  defined  by 


(2.78) 


the  random  errors  for  four  important  parameter  values  may  be 
approximated  as  follows: 


d1  ’xy 
2 


Dv«Q 


(f)| 
h 


/^|Yxy(f)| 


(2.79) 


(2.80) 


c[|  A(f )  0  * 


^-DT^ylfl  | 

^V:IT|qxy(f)l 


(2.81) 


(2.82) 


These  formulas  can  also  be  used  to  determine  the  required  values 
for  n^  to  achieve  acceptable  random  errors  for  these  four  parameter 
values  under  assumed  conditions  for  the  coherence  functions. 

This  concludes  Section  2. 
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3.  NONLINEAR  DRIFT  FORCE  MODELS 


Section  2  details  a  mathematical  procedure  to  determine 
the  spectral  decomposition  of  nonlinear  drift  forces  due 
to  random  ocean  waves  hitting  a  floating  structure  such  as 
a  ship.  More  generally,  it  shows  how  to  determine  the 
spectral  decomposition  of  random  data  passing  through 
parallel  linear  and  nonlinear  systems,  where  the  nonlinear 
system  involves  a  square-law  envelope  detector.  Different 
linear  operations  that  may  be  present  in  the  parallel 
linear  and  nonlinear  paths  can  be  identified  based  solely 
upor.  measurements  of  the  input  data  and  total  output  data. 
The  input  is  assumed  to  be  a  Gaussian  stationary  random 
process  with  arbitrary  autospectral  density  function. 
Formulas  are  stated  to  evaluate  statistical  errors  in 
various  estimates  made  from  measured  data.  These  formulas 
can  also  be  used  to  design  experiments  that  will  achieve 
desired  statistical  errors. 


-59- 


3.1  FORMULATION  OF  DRIFT  FORCE  MODEL 


The  general  problem  to  be  analyzed  for  slowly  varying  drift 
forces  on  floating  structures  is  shown  in  Figure  3.1.  A  ship  is 
assumed  to  be  subjected  to  random  ocean  waves  represented  by 
the  wave  elevation  input  x(t) .  This  produces  a  force  F(t)  on 
the  ship  that  results  in  the  ship  motion  output  y(t) .  The 
spectral  decomposition  of  y(t).  is  desired  that  is  due  to  x(t). 
Random  data  analysis  techniques  from  References  [3,M  will  be 
applied  to  obtain  desired  results. 


Figure  3.1  Illustration  of  drift  force  problem. 

The  force  F(t)  acting  on  the  ship  in  Figure  3.1  is  assumed 
to  consist  of  two  components: 

(1)  a  linear  term  that  is  proportional  to  x(t) , 

(2)  a  nonlinear  term  that  is  proportional  to  the  squared  enve¬ 
lope  signal  of  x(t). 
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This  nonlinear  term  is  called  the  slowly  varying  drift  force, 
References  [1A*1^].  Thus  the  force  F(t)  acting  on  the  ship  can 
be  expressed  as 


F  ( t)  =  k^Ct)  +  k2u(t) 

(3.1) 

x(t)  = 

wave  elevation  input 

(3.2) 

u  ( t )  = 

squared  envelope  signal  of  x(t) 

(3.3) 

l'k2 

proportionality  constants 

(3.1*) 

Equation  3.1  can  be  made  to  apply  to  more  general  situations  by 
letting  k^  and  k2  be  functions  of  frequency.  This  will  be 
done  in  the  following  development. 

The  ship  motion  output  y(t)  in  Figure  3.1  will  be  represent¬ 
ed  by  the  nonlinear  drift  force  model  in  Figure  3.2  such  that  the 
total  output  record 


y  (t) 

=  y1(t)  +  y2(t)  +  n(t) 

(3.5) 

where 

yl(t)  = 

linear  output 

due  to  x(t) 

(3.6) 

y2<t)  = 

uncorrelated 

nonlinear  output  due  to  u(t) 

(3.T) 

n  ( t)  = 

uncorrelated 

zero  mean  Gaussian  output  noise 

(3.8) 

In  Figure  3.2, the  quantities  H^(f)  and  H2(f)  are  frequency 
response  functions  of  constant  parameter  linear  systems  that 
incorporate  the  constants  k^.  and  k2  from  Equation  3.1. 
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x(t) 


y  (t) 


Figure  3.2  Nonlinear  drift  force  model. 


Twc  problems  will  be  treated  using  the  model  of  Figure  3.2. 

(1)  Spectral  Decomposition  Problem.  Given  ( f )  and  H2(f) 
plus  measurement  of  x(t)  only,  determine  the  spectral  pro¬ 
perties  of  y^(t)  and  y2(t).  If  y(t)  is  measured  as 
well  as  x(t) ,  determine  also  the  spectral  properties  of 

n  (t)  . 

(2)  System  Identification  Problem.  From  measurements  of 

x(t)  and  y(t) ,  identify  the  optimum  frequency  response 
functions  H^(f)  and  H 2 ( f )  to  minimize  the  autospectrum 
of  n (t) . 

A  similar  analysis  was  previously  carried  out  in  Section  2  for  the  nonlinear 
wave  force  model  of  Figure  2.^  that. represents  Figure '2.2.  Not  hat  Figure 
3.2  is  a  special  case  of  Figure  1.1*  and  can  be  replaced  by  an  eq  alent 
tvo-input/single-output  linear  model  using  the  inputs  x(t)  and  u(t). 
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3.2  BASIC  FORMULAS  FOR  NONLINEAR  DRIFT  FORCE  MODEL 

In  the  nonlinear  model  of  Figure  3.2,  the  linear  output  term 
y-^t)  is  due  to  x(t)  passing  through  a  constant  parameter 
linear  system  defined  by  the  frequency  response  function  H^(f)  . 
The  function  H^f)  is  the  Fourier  transform  of  a  linear  weight¬ 
ing  function  h^(x)  such  that  in  the  time  domain 


y,(t)  =  f  h1(x)x(t  -  t)  di 

±  J— 00 

(3.9) 

and  in 

the  frequency  domain 

Y]_(f)  =  H1(f)X(f) 

(3.10) 

where 

X(f)  and  Y1(f)  are  Fourier  transforms  of  x(t) 

and 

y1(t) , 

respectively.  Theoretically, 

r  00 

Hl(f)  =  J  h1(i)e"j2lTfTdT 
—00 

/•  0 0  , 

(3.11) 

X(f)  =  J  x  (t)  e"^2irftdt 

J  —00 

(3.12) 

Y.(F)  =  f  y1(t)e'^2ljftdt 

J—oo 

(3.13) 

It  is  assumed  here  and  in  following  formulas  that  mean  values  are 
always  removed  prior  to  computing  Fourier  transforms. 

Ihe  nonlinear  output  term  in  FiSure  3.2.  is  due  to  x(t) 

passing  through  a  square-law  envelope  detector  to  produce  u(t), 
followed  by  u(t)  passing  through  a  constant  parameter  linear 
system  described  by  a  frequency  response  function  ^(f)  .  The 
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function  H 2 (f)  is  the  Fourier  transform  of  a  linear  weighting 
function  h2(x)  where,  in  general,  h2(x)  ^  h^(x)  and 
H2(f)  ^  H^(f).  In  the  time  domain 


y2(t) 


f  h2(i)u(t  -  i)dx 


(3. 1^) 


and  in  the  frequency  domain 


Y2(f)  =  H2(f)U(f)  (3.15) 

where  U(f)  and  Y2(f)  are  Fourier  transforms  of  u(t)  and 
y2(t),  respectively.  Theoretically, 


H2(f)  = 

*  j 

r 

’  —00 

h2(T)e-32,r£TdT 

(3.16) 

U(f)  - 

r  co 

'  —00 

u(t)e'22,,dtat 

(3. IT) 

Y2(f)  = 

r 

/  —00 

, ,  .  -j  2irft . . 
y2(t)e  J  dt 

(3.18) 

As  proved  in  Reference 

[3] 

,  the  output  u(t) 

of  the  square- 

law  envelope  detector  is  given  by 

u(t)  =  x^(t)  +  x^(t) 


where 


x(t)  =  Hilbert  transform  of  x(t)  (3.20) 

When  x(t)  has  a  zero  mean  value,  then  x(t)  will  also  have  a 
zero  mean  value.  However,  the  mean  value  of  u(t)  ,  denoted  by 
E[u(t)],  is 
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(3.21) 


Uu  =  E[u(t)  ]  =  2 ax 


from  Equation  3.19  and  the  fact  that 


E [x2 ( t) ]  =  E [x2 ( t) ] 


(3.22) 


As  usual,  the  notation  E[  J  denotes  the  expected  value  of  the 
quantity  inside  the  brackets.  The  Fourier  transform  of  x(t) , 
denoted  by  X(f)  ,  is 


X(f)  =  f  x(t)e_j2Trftdt 

J  ^  CO 


(3.23) 


How,  Fourier  transforms  of  both  sides  of  Equation  3.19  yield 

/00  ^ 

[X(a)X(f  -  a)  +  X(a)X(f  -  a)] da  (3.21+) 

>00 


Reference  [3]  also  proves  that 


X(f)  =  B  ( f )  X ( f ) 


(3.25) 


where 


B ( f )  =  -jsgnf  = 


f  >  0 
f  =  0 
f  <  0 


(3.26) 


Hence 


U(f)  =  J 


[1  +  B(ct)B(f  -  a)]X(a)X(f  -  a)da 


(3.27) 


and  Equation  3.15  becomes 
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Y2(f)  =  H2(f)  /  u  +  B(a)B(f  -  a)  ]  X(a)  X(f  -  a)dct  C 

Equations  3.10  and  3.28  show  how  to  compute  Y^(f)  and  Y2(f) 
from  knowledge  of  X(f) ,  H^(f)  and  H2(f).  Equations  3. 2k  to 

3.27  show  how  to  compute  U(f)  and  U(f)  from  knowledge  of 
X(f)  . 

The  factor  B(a)B(f  -  a)  in  Equations  3.27  and  3.28  can  be 
determined  as  follows.  From  Equation  3.26,  the  quantity 


~3 

B  (a)  =  0 


a  >  0 
a  =  0 
a  <  0 


Ar  appropriate  plot  for  B(a)  is 


From  Equation  3.2 6,  the  quantity 

-j 

a 

<  f 

B ( f  -  a)  = 

0 

a 

=  f 

■  j 

a 

>  f 

*1 « 


Hence,  for  any  f  >  0 ,  an  appropriate  plot  for  B(f  -  a)  is 


B  (  f  -  a) 


while  for  any  f  <  0,  an  appropriate  plot  is 


B  ( f  -  a) 


Thus,  for  any  f  >  0,  the  product 


B  (a)  B  ( f  -  a)  = 


1  a  <  0 

-1  0  <  a  <  f 

1  a  >  f 


(3.31 


while  for  any  f  <  0,  the  product 
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B  ( a )  B  ( f  -  a ) 


1 

-1 

1 


a  <  f 
f  <  a  <  0 
a  >  0 


(3.32) 


At  values  of  a  -  0  and  a  =  f,  note  that 

B(ot)B(f  -  a)  =  B  ( 0)  B  ( f )  =  0  (3.33) 

Before  applying  these  formulas,  a  brief  discussion  will  be  given 
on  previous  drift  force  models. 


3.3  PREVIOUS  DRIFT  FORCE  MODELS 


Previous  drift  force  models,  such  as  those  assumed  in 
References  [12-143  are  considerably  more  complicated  be¬ 

cause  of  employing  a  bilinear  (quadratic)  weighting  function 
^2  (^l'^2^  and  bilinear  (quadratic)  frequency  response 
H2(f/g)  function  to  represent  the  nonlinear  path.  Specifically, 
Figure3.2  is  extended  to  a  more  general  Figure  3.3  where  y^(t)  is 
the  same  but  y 2(t)  is  replaced  in  the  time  domain  by 

y2(t)  =  J  h2(T1,T2)x(t  -  T1)x(t  -  t  2 )  dx  -^dx  2  (3.3*0 

and  in  the  frequency  domain  by 


Y2(f) 


(a,  f 


a)X(a)X(f  -  a)da 


(3-35) 


Here,  X(f)  and  Y2(f)  are  Fourier  transforms  of  x(t)  and 
y2(t),  respectively,  while  the  bilinear  frequency  response 
function  H2(f,g)  is  the  double  Fourier  transform  of  the  bilinear 
weighting  function  h2(Y^,Y2),  namely, 

/•  -  j  2tt  (  f  x  +  3  x  ) 

h2(x1,r2)e  1  c  dTldT2  (3.36) 


This  bilinear  frequency  response  function  H2(f,g)  in  two 
frequency  variables  f  and  g  is  much  more  difficult  to  compute 
and  interpret  than  the  alternative  H 2 ( f )  in  Figure  3.2.  Mathematical 
developments  of  bilinear  functions  are  contained  in  References  [1,2 ].  Some  books 
that  discuss  these  matters  are  References  [1*4,15,16]. 
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x(t) 


y  ( t) 


Figure  3.3  Parallel  linear  and  bilinear  systems. 


Comparison  of  Equation  3.35  for  Y2  (f)  in  Figure  3.3  with 
Equation  3.28  for  Y2(f)  in  Figure  3.2  shows  that  the  same  Y2(f) 
will  occur  if 

H2(a,f  -  a)  =  H2(f)[l  +  B(cOB(f  -  a)]  (3.37) 

Thus,  knowledge  of  the  linear  frequency  response  function  H2(f) 
in  Figure  3.2  can  yield  an  equivalent  bilinear  frequency  response 
H2(a,f  -  a)  for  Figure  3.3.  However,  knowledge  of  the  bilinear 
frequency  response  function  H2(a,f  -  a)  will  not  yield  an 
equivalent  linear  frequency  response  function  H2(f)  for  Figure  3.2 
To  explain  this  matter,  note  from  Equations  3.31  and  3.32 
that  the  function  [1  +  B(ct)B(f  -  a)]  will  be  zero  for 
0  <  a  <  f  when  f  >  0,  and  will  be  zero  for  f  <  a  <  0  when 
f  <  0.  For  these  values  of  a  and  f,  one  cannot  solve 
Equation  3.37  for  H2(f)  to  satisfy  an  arbitrary 
H2(a,f  -  a)  ^  0.  Even  when  the  factor  fl  +  B(a)B(f  -  a)}  f  0, 
from  Equation  3.37,  H2(f)  must  satisfy 
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H2 (a, f  -  a) 


(3.38) 


H  (f)  =  - 

[1  +  B(a)B(£  -  a)] 

This  requires  the  right-hand  side  to  be  independent  of  flf ,  a 
situation  that  is  highly  unlikely  to  occur  in  practice  for  an 
arbitrary  H2(a,f  -  a)  .  The  conclusion  is  that- one  cannot  use 
past  results  about  H2(a,f  -  a)  to  determine  an  equivalent 
H2(f)  for  Figure  3.2.  Instead,  future  work  should  be  based  on 
new  results  obtained  by  following  the  procedures  described  in 
this  report. 
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3.-  SPECTRAL  DECOMPOSITION  PROBLEM 

For  Figure  3.2, the  Fourier  transform  of  Equation  3.5  proves 

that 

Y  ( f )  =  YL(f)  +  Y2(f)  +  N(f)  (3.39) 

where 


Y1(f) 

=  H1{f)X(f) 

(3.b0) 

Y2(f) 

=  H2(f)U(f) 

(3. hi) 

N  ( f ) 

=  Fourier  transform  of  n(t) 

(3.h2) 

Equation  3.39  will  be  used  henceforth  in  place  of  Equation  3.5 
and  all  further  analysis  will  be  done  with  frequency  domain 
quantities . 

As  derived  and  illustrated  in  References  [3,1*],  from  knowledge 

of  X(f)  and  'Y(f)  for  a  number  of  different  stationary  random  records  of 

length  T  ,  their  one-sided  (f  >  0)  autospectral  and  cross- 

spectral  density  functions  can  be  computed  by  the  direct  form¬ 
ulas 

GxX(f)  =  |  E[X*(f)X(f)  ]  (3.1*3) 

Gyy(f)  =  I  E[Y*(f)Y(f)l  (2.kh) 

Gxy(f)  c  |  E[X*(f)Y(f)  ]  (3.«*5) 

This  procedure  does  not  require  the  computation  of  any  associated 
autocorrelation  or  cross-correlation  functions. 
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From  Equations  3.39  to  3.1*2,  the  total  output  autospectral 
density  function  is  given  by  the  formula 


Gyy(f)  *  Vi"’  +  V°  +  Gnn<f) 

(3.1*6) 

where 

Vi1"’ =  |Hi<f)|2v(f) 

(3.1*7) 

Gy2y2<f)  =  |H2(f)|2GUu(f) 

( 3 .  h  R ) 

G  (f)  =  autospectrum  of  n(t) 

(3A9) 

The  cross-spectral  density  function  between  y ^  ( t )  and 

y2(t)  , 

namely. 

Gy1y2(f)  =  HI(f)H2(f)Gxu(f>  -  0 

(3.50) 

because 

x(t)  being  Gaussian  makes 

Gxu(£)  -  0 

(3.51) 

Also,  assumptions  about  n(t)  make 


V(fl  =  =  0  tJ'52) 

The  spectral  quantity  'G  u(f)  in  Equation:  3.  1*8  can  be  computed 
knowledge  of  U(f)  for  a  number  of  different  stationary  random 
records  of  length  T  by  the  direct  formula 

Guu(f)  =  §  E[U*(£)U(«]  (3.53) 

Since  U(f)  is  known  from  X(f)  ,  Equation  3.53  actually  shows  that 
Guu(£)  is  a  function  of  X(f) . 


-73- 


A  useful  theoretical  formula  is  derived  in  Section  3.8.1  to 
compute  Guu(f)  from  knowledge  of  the  input  autospectral  density 
function  G  (f)  .  For  any  f  >  0,  the  formula  is 

AA 


Equation  3.51*  is  applicable  to  both  the  spectral  decomposition 
problem  and  the  system  identification  problem,  while  Equation  3-56 
is  applicable  only  to  the  system  identification  problem. 

The  spectral  decomposition  problem  has  now  been  solved. 

From  measurement  only  of  x(t),  one  can  compute  G  (f)  and 

aA 

Guu(f).  Then  from  knowledge  of  H^(f)  and  ^(f),  one  can 

compute  the  linear  part  of  the  output  autospectrum  G  (f)  due 

ylyl 

to  the  wave  elevation  force  by  Equation  3.**T  ,  and  the  nonlinear 

part  of  the  output  autospectrum  G  (f)  due  to  the  nonlinear 

y  2y2 

drift  force  by  Equation  3.^8  . 


If  y(t)  is  measured  as  well  as  x(t) ,  then  one  can 


compute  Gyy ( f ) .  It  is  now  possible  to  obtain  the  output  noise 
autospectrum  Gnn(f)  from  Equation  3.^6  by  the  formula 


-  Gyy(f)  -  Vl'*’  ‘  GW£) 


(3.57) 


Note  also  that  if  Gnn(f)  is  zero,  then  G  (f)  can  be  Predicted 
without  measuring  y(t)  by  the  formula 


Gyy(f)  -  <*>  ♦  V2(f) 


(3.58) 


Coherence  functions  are  required  to  evaluate  statistical 
errors  in  estimates  to  be  given  in  Section  3.7.  The  linear  coher¬ 
ence  function  y *  (f)  between  the  input  x(t)  and  the  total 

xy 

output  y (t)  is  defined  by 


Y*(f>  = 


lGxy(f>T 

GXX(f>Gyy(f) 


Gyiyi(f) 
G  (f) 

yy 


(3.59) 


This  gives  the  percentage  of  the  total  output  autospectrum  due 


to  the  linear  operations  on  x(t).  A  nonlinear  coherence  function 
2 

q  (f)  can  be  defined  here  by 
xy 


(f) 


(3.60) 


This  gives  the  percentage  of  the  total  output  autospectrum  due 
to  the  nonlinear  operations  on  x(t).  In  terms  of  these  two 
coherence  functions,  the  output  noise  autospectrum  of  Equation 
3-57  becomes 
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(3.61) 


G  (f)  =  fl  -  Y2  (f)  -  q  2  (f)l  G  (f) 
nn  L  '  xy  xy  J  yy 

Figure  3.2  will  be  a  good  nonlinear  drift  force  model  at  those 
frequencies  where  Gnn(f)  is  close  to  zero,  namely,  where  the 
2  2 

sum  [y  (f)  +  cv„(f) ]  is  close  to  unity.  Otherwise,  Figure  3.2 
xy  xy 

is  a  poor  model. 
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3-5  SYSTEM  IDENTIFICATION  PROBLEM 

Assume  now  that  the  system  properties  of  H^(f)  and  H2(f) 
are  not  known  in  Figure  3.2.  Optimum  properties  are  to  be  identified 
from  measurements  of  x(t)  and  y(t)  based  upon  minimizing  the 
autospectrum  G  (f)  °f  n(t)  over  all  possible  choices  of 
linear  systems  to  predict  y(t)  from  x(t) .  It  is  assumed  as 
before  that  x(t)  follows  a  Gaussian  distribution  with  zero 
mean  value.  Nonzero  mean  value  inputs  are  treated  in  Section  3.6. 

In  Figure  3.2, from  Equations  3.39  to  3.12,  the  basic  Fourier 
transform  relation  is 

Y(f)  =  H]_(f)X(f)  +  H2(f)U(f)  +  N(f)  (3.62) 

where  X(f),  U(f)  and  Y(f)  can  be  calculated  from  the  given 

x(t)  and  y(t).  From  Equation  3.62,  the  cross-spectral  density 

function  G  (f)  between  x(t)  and  y(t)  satisfies 
xy 

Gxy(£>  '  Hl(£)Gxx(f)  13  63) 

provided  that 


Gxu(f> 


0 


G _ (f)  =  0 


(3.61) 

(3.65) 


Equation  3.61*  occurs  because  x(t)  and  u(t)  will  be  uncorre¬ 
lated  when  x(t)  is  Gaussian  data.  Equation  3.65  occurs  by 
assuming  x(t)  and  n(t)  are  uncorrelated. 

From  References  [3,1*],  the  optimum  linear  system  H  (f) 
between  x(t)  and  y(t)  is  given  by  the  formula 


x  (t)  = 

0  Gxx<£) 


(3. 66) 


For  Figure  3.2, from  Equation  3.63, 


Gxy(f)  ,  ,  -j<jJl(f) 

H,  (f)  =  -  =  |  H.  (f)  j  e 


Gxx<f> 


(3.67) 


Thus  H1(f)  is  the  same  here  as  the  optimum  linear  system 

H  (f)  that  produces  the  minimum  at  all  f  over  all 

o  nn 

possible  linear  systems.  The  optimum  system  automatically  makes 
n(t)  uncorrelated  with  x(t)  so  that  Equation  3.65  applies 
without  the  necessity  to  assume  in  advance  that  x(t)  and  n(t) 
will  be  uncorrelated.  Equation  3.67  shows  that  the  complete  H-^f) 
in  both  gain  and  phase  can  be  identified  using  only  x(t)  and 
y (t)  . 

From  Equation  3.62,  the  cross-spectral  density  function 
GUy(f)  between  u(t)  and  y(t)  satisfies 


Guy<f>  *  H2(f)Guu(£) 


(3.60) 


provided 


G  !  =0 

ux 


(3.69) 


Gun<f>  =  0 


(3.70) 
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These  two  results  follow  from  Equations  3. 61*  and  3.65.  They  also 
occur  because  x(t)  and  u(t)  will  be  uncorrelated  when  x(t) 
is  Gaussian  data,  and  because  n(t)  and  u(t)  will  be  automati¬ 
cally  uncorrelated  when  H^(f)  satisfies  Equation  3.67.  For 
Figure  3. 2,  from  Equation  3.68  , 


Guy(f)  ,  , 

H  (f)  =  -22 -  =  |H,(f)  |  e 

G  (f) 
uu 


(3.71) 


In  practice,  G  (f)  can  be  computed  directly  by  Equation  3*53 


and  G  (f)  can  be  computed  directly  by 


Guy(f)  =  I  E[U* (f) Y (f ) ] 


(3.72) 


Since  u(t)  is  known  from  x(t),  Equations  3.53,  3.71  and  3.72  show 

that  the  complete  H2(f)  in  both  gain  and  phase  can  be  identified 

using  only  x(t)  and  y(t).  Alternate  theoretical  ways  to 

compute  Guu(f)  and  Gu^(f)  are  9iven  in  Equations  3.5^  and  3.56, 

bwsed  upon  derivations  in  Sections  3-8.1  and  3.8.2.  A  different  new 
theoretical  way  to  identify  H2(f)  is  derived  in  Section  3.8.3 

After  H^(f)  and  H2(f)  have  been  computed  by  Equations  3.67 

and  3.71  ,  respectively,  the  spectral  decomposition  problem  can  then 

be  solved  using  Equations  3.1*7  ,  3.1*8  ,  and  3.57.  To  evaluate 


statistical  errors  in  estimates,  one  should  also  compute  the  linear 


coherence  function 
herence  function  q 


2 

Y 

'xy 


(f) 


by  Equation  3.59  and  the  nonlinear  co- 
by  Equation  3.60. 
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3.6  NONZERO  MEAN  VALUE  INPUT 

Suppose  that  the  Gaussian  input  x(t)  has  a  nonzero  mean  value  given  by 

u*  =  £[x(t)}  f  0  (3.73) 

Then  x(.t)  can  be  expressed  as 

x(t)  =  lx(t)  -  yx]  +  yx  (3.7M 

consisting  of  a  variable  input  term  [x(t)  -  p  ]  with  zero  mean  value,  plus 
a  constant  term  u  .  Here  x(t) ,  the  Hilbert  transform  of  x(t) ,  will 

A 

still  have  a  zero  mean  value  since 

-  E[x(t)]  =  B(0)px  =  0  (3.75) 

and  B(0)  =  0  from  Equation  3.26.  In  place  of  Equation  3.19,  the  output  of 
the  square- law  envelope  detector  Is  now  given  by  w(t)  where  w(t)  is 


2  2 

W(t)  =  { [  X  (t)  -  yj  +  Px}  +  X  (t) 

=  [x(t)  -  Ux]2  +  X2(t)  +  2px[x(t)  -  yj  +  Ux 
=  u ( t )  +  2ux[x(t)  -  px]  +  mx  (3.76) 
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The  last  term  jU^  produces  an  effect  in  the  spectrum  of  w(t)  only  at 
f  =  0.  For  £  f  0,  Equation  3.76  consists  of  a  nonlinear  square-law  enve¬ 
lope  detector  in  parallel  with  a  linear  operation  2y  as  shown  in  Figure  3.4. 

A 

Note  that  Figure  3.*+  simplifies  to  the  square- law  envelope  detector  by  itself 
when  ux  =  0. 


Figure  3.4  Nonlinear  system  when  yx  f  0. 

Return  now  to  the  original  problem  covered  by  Figure  3.2.  For  situations 
where  the  input  changes  to  [x(t)  -  y  ]  with  y  f  0  ,  Figure  3.2  becomes 
the  new  extended  Figure  3-5.  In  Figure  3.5,  the  two  outputs  y^(t)  and  y3(t) 
will  be  correlated  because  of  the  linear  operation  in  the  nonlinear  path. 
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[x(t) 


3.5  Nonlinear  drift  force  model  when  ux  f  0. 


The  optimum  linear  system  HQ(f)  in  Figure  3.5  is  computed  as  before  by 
the  equation 


H  (f) 
o 


yfl 

yo 


Hj(f)  *  2:jxH2Cf) 


(3. 


Note  that  H  (f)  is  now  a  function  of  both  H.  (f)  and  Figure  3.5 

should  now  be  revised  so  that  H^(f)  is  replaced  by  H  (f)  to  produce  a 
new  output  y4(t)  in  place  of  y^(t).  The  output  y3(t)  in  Figure  3.5  will 
then  become  the  previous  output  y2(t)  in  Figure  3.2  and  will  preserve  the 
sum 


yAt)  +  y3(t)  =  y4(t)  +  y2(0 


(3. 


This  will  give  Figure  3.6  where  the  two  outputs  y^(t)  and  y^Ct)  will  be 

uncorrelated  because  of  the  properties  of  H  (f)  and  the  Gaussian  nature  of 

o 

the  input  [x(t)  -  p  ] . 


y(tj 


Figure  3.6  Equivalent  model  to  Figure  3.5. 

The  following  Fourier  transform  relations  are  satisfied  by  the  quantities 
shown  in  Figure  3.5  and  3.6.  Various  applications  are  feasible  depending  upon 
what  is  known  and  what  can  be  measured. 


Also  ,  for  f  ^  0, 


X(f) 

*  2T[xCt) 

‘  vx] 

(3.79) 

U(f) 

*  3f[u(t) 

'  »u] 

(3.80) 

W(f) 

■  ?[w(t) 

(3.81) 

YCf) 

=  STlyCt) 

■  V 

(3.82) 

WCf)  ’ 

=  U(f)  «• 

2yx  X(f) 

(3.83 
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(3.8*0 


Y1Cf)  -  JtYiCt)  -  u  1  =  H^nxcf) 

Y2(f)  =  ?[y2(f)  -  ji  ]  -  H2(f)U(f)  (3.85) 

Y3Cf)  -  7[yz(t)  -  Uy^}  =  H2(f)W(f)  (3.86) 

Y4C£)  -  ?[y4Ct)  -  U  ]  =  H0Cf)X(£)  (3.87) 

where  H  (f)  is  given  by  Equation  3.77.  The  output  term 

Y(f)  *  Y4(f)  +  Y2Cf)  +  N(£)  (3.88) 


From  Equation  3.88,  the  output  autospectral  density  function  is 


V*  -  wf)  *  v2(f)  *  G-(f)  <3-69> 

where 

V«C£)  ’  'H0Cf)'2G>«C£)  <3-90) 

Gy  (0  -  |H2(f)|2Guu(f)  (3.91) 

G  (f)  =  autospectrum  of  unmeasured 
!  extraneous  output  noise 


From  knowledge  of  H^(f) ,  H2(f)  and  x(t) ,  one  can  compute  Hc(f),  G^tf) 
and  G^Cf)  so  as  to  be  able  to  predict  Yy  y  (f)  and  Gy  y  (f).  When 
G^tf)  =  0,  this  will  predict  the  total  output  spectrum  Gyy(f). 

To  identify  the  optimum  system  properties  of  HQ(f),  H^(f)  and  H2(f) 
from  measured  input/output  data,  the  appropriate  equations  to  use  are 


*  Ho<«Gxx<« 

(3.92) 

WCf)  ■  ¥0Guu(« 

(3.93) 
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As  in  Section  3.5,  from  measurements  of  x(t)  and  y(t),  one  can 
compute  the  spectral  quantities  shown  in  Equations  3.92  and 3.93  .  This 
gives 


,  cm 

H0(f)  =  - 

(3.9b) 

H2(f)  = 

Guuf« 

(3.95) 

Then,  from  Equation  3.77 

Hx(f)  =  H0Cf)  -  2yxH2Cf) 

(3.96) 

Other  spectral  relations  of  possible  interest  for  special  applications 

when  y  f  0  are 

A. 

Vi(f)  =  |Hit£)|2G**f£) 

(3.97) 

■  IH2Cf)|2(^Cf) 

(3.98) 

’  Guutf)  ‘  Gxx<« 

(3.99) 

Vj'O  " 

(3.100) 

V/*  ’  ° 

(3.101) 

This  completes  the  analysis  for  nonzero  mean  value  inputs. 
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3.7  STATISTICAL  ERRORS  IN  ESTIMATES 


Return  back  to  Figure3.2and  the  equations  in  Sections  3.1* 
and  3.5.  For  this  error  analysis,  replace' Figure  3.2  by  a  two-input/ 
single-output  linear  model  with  uncorrelated  inputs  x(t)  and  u(t)  as  shown 
in  Figure  3.7. 


x  ( t ) 


u(t) 


Figure  3.7  Two-input/single-output  linear  model 
for  Figure  3.2. 

Statistical  random  error  formulas  will  now  be  stated  for 
estimates  of  the  following  six  quantities.  These  formulas  are 

derived  in  Reference  [3]  and  summarized  in  Reference  [  U] . 


(1) 

8v  y  (£) 
yiYl 

computed 

by 

Equation 

3.1*7 

(2) 

^y2y2(f) 

computed 

by 

Equation 

3.1*8 

(3) 

a2  , 

Yxy(f) 

computed 

by 

Equation 

3.59 

(4) 

?xy(f> 

computed 

by 

Equation 

3.60 

(5) 

(f ) 

computed 

by 

Equation 

3.67 

(6) 

H  2  {  f ) 

computed 

by 

Equation 

3.71 
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The  normalized  random  errors  for  the  two  output  spectra 

A  A 

estimates  G  (f)  and  G  (f)  are  approximated  by 
ylyl  y2y2 


[2  -  y2  ( f ) ] 1/2 

elG„„(f)]«  - ^ - 


ylyl 


(3.102) 


I9xy(f)l^ 


E[Gy2y2(£)1=! 


[2  -  qxy (£) ] 1/2 


(3.103) 


where  n^  is  the  number  of  independent  averages  used  to  compute 
the  original  spectral  density  quantities.  Required  values  of 


n,  to  achieve  desired  values  of  e[G  ( f ) ]  or  e[G  (f)] 

ylyl  y2y2 


are  obtained  by  solving  Equations  (73)  and  (74)  for  nd> 

The  normalized  random  errors  for  the  two  coherence  function 
'2  *•  2 

estimates  YxvCf)  and  cv„(f)  are  approximated  by 

xy 


2  ^  t1  -  V(f)l 

e[>xy(f)]»  - ^ - 

l9Xy(f) I  /5d 


(3.101*) 


*2  ^  t1  "  xy ( f } ^ 

e[c^y(f)]»  - & - 


|qxy(f)| 


Bias  errors  for  these  estimates  are 


(3.105) 


~2  [1  "  yxv(f)]2 

b[Yxy(f)l  »  - ^ - 


(3.106) 
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T [4^  (f )  ]  ~  sin-1  |  e  [  |  H1  ( f )  |  ]  |  (3.110) 

^■[^(f)]  «  sin'1  |e[fH2(f)f]|  (3.111) 

For  small  values  of  e(|H^|]  and  e[|H2|]  ,  Equations  3.110  and 
3.111  can  be  further  simplified  to 
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(3.112) 


of^Cf)]  &  £[|H1(f)  |  ] 
a[$2(f)  ]  C  e[  |  H2  ( f )  |  ]  (3.113) 

Various  bias  error  in  frequency  response  function  estimates  are 
discussed  also  in  Reference  [3]  that  should  be  minimized  as 
much  as  possible  in  practice. 

3.8  DERIVATIONS  OF  THEORETICAL  FORMULAS 

To  make  this  report  more  self-contained,  some  special  theoretical 

formulas  are  derived  here  that  supplement  material  discussed  in  earlier  sections. 

These  formulas  deal  with  quantities  G  (f),  G  (f)  and  H  (f)  that  are  involved 

uu  uy  2 

in  the  nonlinear  drift  force  models  of  Figures  3.2  and  3.6.  Readers  not  interested 
in  these  matters  should  proceed  to  Section  4. 
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3.8.1  Theoretical  Formula  for  G  (f) 

uu 


A  useful  theoretical  formula  will  now  be  derived  to 
determine  G  ,(f)  from  knowledge  of  G  (f).  From  Equation  3.2l, 

Till  XX 


0(f)  =  r 

J  —  c 


[X(cOX(f  -  a)  +  X(a)X(f  -  a)  ]  da 


tl*(f)  =  f  [X*(8)X*(f  -  6)  +  X*  c  6)  x*  Cf  -  6)  ]de 

J  -  oo 


(3. ni*) 


(3.115) 


The  two-sided  autospectral  density  function  S  (f)  of  u(t) 
is  defined  by 


i  E[0*(f)UCf)l 


CB)X*(f-6)  +  X*(8)X*(f-6)}{X(a)X(f-a) 
+  X(p)X(f  -  Q  )}  ]dad|3 


(3.116) 


In  Equation  3.116  the  integrand  involves  a  total  of  four  different 
fourth-order  moments  of  Gaussian  data  that  breaks  down  into 
products  of  pairs  of  second-order  moments.  A  typical  fourth- 
order  moment  can  be  replaced  by  the  three  product  pairs 


E[X*(B)X*(f  -  6)X(a)X(f  -  a)]  =  E[X*(B)X*(f  -  8)  ]E[X(a)X(f  -  a)  ] 


(3-117 ) 


+  E[X*(B)X(a)  ]E[X*(f  -  8)  X(f  -  a)  ]  +  E[X*(8)X(f  -  a)]E[X*(f  -  6)X(a)] 


The  first  pair  of  second-order  moments  in  Equation  3. 117  are  given  by 

E[X* (8) X* (f  -  8)]  =  Sxx(B)<51(f)  (3.118) 

E[X  (a) X  (f  -  a)]  =  Sxx(a)61(f)  (3.119) 

where  S  (a)  and  S  (8)  are  two-sided  autospectral  density  functions 

2QC  XX 

of  x(t)  and  where  6^(f)  is  the  finite  delta  function  defined 
by 


61(f) 


61(-f) 


=  T 


(-1/2T)  <  f  <  (1/2T) 


(3.120) 


=  0 


otherwise 


&1(f)df  = 


pl/2T) 

(-1/2T) 


61(f)df  =  1 


(3.121) 


The  first  pair  gives  a  contribution  to  Suu(f)  of  Equation  3.116 
only  at  f  =  0  since 


1 

T 


Sxx(a) sxx ( 8) (f ) dad8  = 


(3.122) 


Here 


2 


Sxx(c.)da 


/”  SxX'6)de 


(3.123) 


The  second  pair  of  second-order  moments  in  Equation  3.117  are 
given  by 


E[X*(8)X(a)  ]  =  Sxx(B)  6x(a  -8) 


(3.12U) 
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(3.125) 


E[X  (f  -  8)  X  ( f  -  a)]  =  Sxx(f  -  8)<S1(a  -  8) 


This  second  pair  gives  a  contribution  to  Suu(f)  of  Equation  3.116 
represented  by 


UJ 


Sxx(8)Sxx(f  -  6) 6*( a  -  8)d6da 


Sxx(a)sxx(f  -  a)da  (3.126) 


The  third-pair  of  second-order  moments  in  Equation  3.117  are 


given  by 


E[X  ( 8)  X(f  -  a)  ]  =  Sxx(B)61(f  -  a  -  6) 


E[X*(f  -  8)  X (a)  ]  =  Sxx(f-(3)  £  (f  -  <*-£) 


(3.127) 


(3.128) 


This  third  pair  gives  a  contribution  to  suu(f)  of  Equation  3.116 


represented  by 


n:i 


Sxx(8)sxx(f  -  8)  (f  -  a  -  8)dBda 


-r 


Sxx(f  "  a)Sxx(a)da  (3-129> 


Equation  3-129  is  exactly  the  same  as  Equation  3.126. 
Equations  3.117  to  3.129  prove  for  f  /  0  that 


i  J  Je[X*  (6>)  X*  (f  -  8)X(a)X(f  -  a)  ]  dad  8  ~  2 


Sxx(a)Sxx(f  "  a)da  ( 3. 130) 


Similarly,  one  can  prove  for  f  ^  0  that  the  other  fourth-order 
moments  in  Equation  3.116  lead  to  the  results 

i  f  •  [E[X*(8)X*<f  -  8)X(a)X(f  -  a)  ]dBdct  =  2  f  Sxx  (a)  S  ( f  -  a)  da  (3.131) 

•  —  00  J  J —CO 

till  E  [X*  ( B)  X*  (  f  -  6)  X  (a)  X  ( f  -  a)  ]  d$da  =  2  J  Sxx  (a)  Sxx  (  f  -  a)  da  (3.132) 
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E[X*  (6)  X*  (f  -  6)X(a)X(f  -  a)  ]dSda  =  2  J  Sxx  (a)  &xx  (  f  -  a)  da  (3.133) 

The  derivation  of  Equations  3.131  to  3.133  requires  formulas 
from  Table  13.2  of  Reference  [2].  In  particular 


S  (cO  =  S„  (a) 

XX  xx 


-S_  (a) 

XX 


(3.131) 


where  S  (a)  is  the  Fourier  transform  of  R  (x),  and  R  (x) 

XX  XX  XX 

is  the  Hilbert  transform  of  the  autocorrelation  function  R  ft)  . 

xx v  J 

Hence 


Sxxta)  =  Bfa)Sxxta) 


(3.135) 


where  B(a)  satisfies  Equation  3.29. 

Substitution  of  Equation  3.117  into  Equation  3.116  followed  by 
substitutions  from  Equations  3.130  to  3.133  now  proves  the  theore¬ 
tical  formulas  that  for  any  f  t  0  , 

/oo 

tS 6xx(“)Sxx(f-a)  +5xx<c,)rxx(f  <5-136) 

>00 

Suu(f)  =4  /  [1  +  B(a)B(f -a)]Sxx(a)S  (f -a)dd  (3.137) 

J  _00  u 


These  formulas  show  how  to  compute  Suu(f] 


Sxx(f>- 


from  knowledge  of 


The  one-sided  autospectral  density  functions  Guu(f)  and 

6  (f)  are  related  to  the  two-sided  autospectral  density  functions 

Suu(f)  and  Sxx(f)  by 


Guu(f’ 


2S  (f) 
uu 


0 


f  >  0 
f  <  0 


(3.138) 
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G  (f)  = 
xx 


2S  (f) 
xx 


f  >  0 


f  <  0 


(3.139) 


Also 


Suu(f)  =  suu(-f)  ~  I  Guu(f^ 


( 3.H0) 


Sxx(f)  -  Sxx(-f)  -  J  Suu(f> 


(3.H1) 


Consider  any  value  of  f  >  0  and  apply  Equation  3.31  to 
Equation  3. 137.  Now 


Suu<£>  ■  8 


/  Sxx<a>sxx“  -“>da+  8  f  Sxx(t“)sxx'£  "  “)da  (3-lh2> 

J  —00  J  P 


This  is  the  same  as 


Suu<f>  *  8  f”  Sxx(-B)Sxx(f  +  8)d8  +  Sf  Sxx(a)Sxx(f  "  “,da  i3'li3) 

by  letting  8  =  -a  and  d8  =  -da  .  Replacing  8  by  a  again 
and  using  the  fact  that  S  (-a)  =  S  (a)  ,  S  (f  -  a)  =  S  (a  -  f)  , 


Suu(f) 


-j: 


0  Sxx(a)Sxx(f +  a)da  +  8Jf  Sxx(a)Sxx(“-f)da 


(3.H1) 


For  any  f  >  0  ,  one  now  has  the  equivalent  formula  in  terms  of 
one-sided  autospectral  density  functions 


Guu(f>  Gxx(“)Gxx(f  +  a)d“  +  4Jf  GxxWGxx'“-£>a“ 


(3. 115) 


Equation  3.ll5  shows  how  to  compute  Guu(f)  from  knowledge  of 
Gxx<f>- 
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3.8.2  Theoretical  Formula  for  G  (f) 

_ uy 

Another  theoretical  formula  of  interest  is  the  cross-spectral 
density  function  between  u(t)  and  y(t).  The  two-sided  cross- 
spectral  density  function  S  (f)  is  defined  by 

Suy(f)  =  ^  E  [U*  ( f )  Y  ( f )  ]  (3.11*6) 

where  U(f)  and  Y(f)  satisfy  Equations  3. 27  and  3.39,  respectively. 
Hence 

C  00 

Suy(f)  =  fj  [1  +  B*(a)B*(f  -  a)  ]E[X*  (a)  X*  (f -a)Y(f)]da  (3.H*7) 

For  all  values  of  a  and  f,  one  can  verify  that  the  product 
quantity 

B*  (a)  B*  ( f  -  a)  =B(a)B(f-a)  (3.ll*8) 

The  two-sided  cross-bispectrum  (also  called  the  second-order 
cross-spectral  density  function)  between  x(t)  and  y(t)  is 
defined  for  all  possible  f  and  g  by  the  two-dimensional  func¬ 
tion 

Sxxy(f'g)  *  k  E[X*(f)X*(g)Y(f  +  g)] 

As  a  special  case,  when  f  =  ot  and  g  =  (f  -  0 1  ),  one  can  compute 

Sxxy(ot'f  _  a)  =  f  E[X*(a)X*(f  -  a)Y(f)  ] 

This  result  will  not  be  zero  for  the  Y(f)  of  Equation  3.39. 
Substitution  of  Equations  3. lU8  and  3.150  into  Equation  3 . 1^7 
yields  the  formula 


(3.11*9) 


(3.150) 
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s,.,.(f)  =  f  [1  +  B (a) B (f  -  a) ] S  (a, f  -  a) da 
U1  J_  00  X  y 


(3.151) 


This  formula  shows  how  to  compute  Suy(f)  from  knowledge  of 
SXxy(a'£-a)- 

The  one-sided  cross-spectral  density  function  G  (f)  is 
related  to  the  twc-sided  cross-spectral  density  function  S  (f)  by 


Guy(f>  " 


2Suy(f) 


f  >  0 
f  <  0 


(3.152) 


Define  the  one-sided  cross-bispectrum  Gxxy  (°*/ f  -0t)  by 


G  ( 
xxy 


/ 

a,  f  -  a)  =  • 


23xxy(a'£-tl)  f  21 

o  f  <  0 


(3.153) 


Now,  for  any  f  >  0  ,  Equation  3.151  becomes 


uy(f)  =  J  [1  +  B(a)B(f  -  a)]  Gxxy(a,f  -  a)  da 


(3.15M 


Applying  Equation  3.31  to  Equation  3.151*  yields 


G(f) 

Uy 


— o 


G _ „(a,f-a)da  +  2  Gvv„  (a  ,  f  -  a)  da 


(3.155) 


Also 


~r 


G^„„(-a,f  +  a)da  +  2  /  G _ (a,f-a)da  (3.156) 


The  cross-bispectrum  functions  in  Equations  3.155  and  3.156  can 


be  determined  for  any  f  >  0  by  the  computations 
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Cxxy(_a,f  +  a)  =  f  E[X*(-a)X*f  +<x)Y(f)  ]  E  [X  fa)  X*  (  f +a)  Y  (  f )  ] 

(3-157) 


Gxxy(a,f-a)  =  |  E[X*(a)X*(f  -  a)Y(f)  ]  E  [X*(a)  X(o- £ )  Y  (£)  ] 

(3.158) 


Equations  3.15^  to  3.158  show  appropriate  ways  to  compute  G  (f) 
from  knowledge  of  G  (a,f  -  a)  . 

Aaj 
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3.8.3  Identification  of 


A  different  new  theoretical  way  to  identify  H2(f)  is 
also  possible  from  measured  x(t)  and  y(t)  by  using 
Equation  3.68  together  with  Equations  3.137  and  3.151.  In 
terms  of  two-sided  spectral  density  functions,  Equation  3.68  is 


Suy(f)  =  H2(f)Suu(f) 


(3.159) 


Substitution  from  Equation  3.137  now  gives 


Suy(f)  =  4H2(f) 


I  [1  +  B  (a)  B  (f  -  a)  ]  Sxx  (a)  Sxx  (f  -  a)  da  (3.l6o) 


However,  Equation  3.151  states  that 


Suy(£’ 


■/* 


[1  +  B(a)B(£  -  a) ]Sxxy(a,f  -  a)da 


(3.161) 


Hence,  one  obtains  the  useful  relation 


W“'f-a)  =  4H2(f>Sxx(o,Sxx(f  “«> 


(3.162) 


This  provides  an  alternate  way  to  estimate  H2(f)  from  computa¬ 
tion  of  the  other  quantities.  Consider  a  change  in  variables 
where  a  =  g  and  f  =  2g.  Then  (f  -  a)  =  g  and  Equation  3.162 
becomes 


Sxxy<9,g) 


-  4H2(2g)S^x(g) 


(3.163) 


?his  is  the  same  as 
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Sxxy(£'f)  =  4H2<2£»s;x(fl 


(3.3.61*) 


Tnus  H2(2f)  satisfies  here  the  simple  relation 


H2(2f) 


xxy 


4Six<£> 


(3.165) 


where 


sxx(f)  =  f  E[X*(f)X(f)  ] 


(3.166) 


SxxyC£»f)  =  T  E[X*Cf)x"(f)Y(2£)] 


(3.167) 


with  Y(2f)  given  by  Equation  .  Finally,  Equation 


proves  that 


«2(f) 


SXXy(f/2,f/2) 


4Sxx(f/'2) 


(3.168) 


This  concludes  Section  3. 
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4.  NONLINEAR  NONSYMMETRICAL  SYSTEMS 


Section  4  develops  various  input/output  relationships  when 
random  data  passes  through  three  types  of  nonlinear 
nonsymmetrical  systems:  (A)  Three-Slope  Systems,  (B) 
Catenary  Systems  and  (C)  Smooth-Limiter  Systems.  Included 
are  formulas  for  input/output  probability  density 
functions,  correlation  functions  and  spectral  density 
functions.  These  particular  nonlinear  nonsymmetrical 
systems  are  of  interest  for  NCEL  applications.  Spectral 
formulas  are  stated  here  using  two-sided  spectral  density 
functions. 
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.1  THREE- SLOPE  SYSTEMS 


A  first  example  of  a  nonlinear  non symmetrical  system  is  a 
three-slope  system  as  sketched  in  Figure  U.l.  Let 


y  =  g(x) 


x 

-A  +  a (x  +  A) 
B  +  b(x  -  B) 


X  £  -A 
x  >  B 


(U.  1) 


For  definiteness,  assume  here  that  the  constants 


0  <  A  <  B ,  0  <  a  <  1, 

The  inverse  relations  are 

x  =  g_1(y)  =  y 

=  -A  +  £(y  +  A)  /a} 

=  B  +  [(y  -  B)/b] 

The  derivative  relations  are 

£ 

=  a 

-  b 


b  >  1  (h.2) 

-A  C  y  £  B 

y  <,  -A  (4.3) 

y  >,  B 

-A  _<  X  £  B 

x  <  -A  (k.k) 

X  >  B 


Thus  the  derivative  is  discontinuous  at  x  =  -A  and  at  x  =  B. 

Let  p(x)  be  the  input  probability  density  function  and 
P2(y)  be  the  output  probability  density  function.  For  this 
problem,  from  Equation  1.5, 

p2(y)  =  H.5) 

| dy/dx | 

where  x  on  the  right-hand  side  should  be  replaced  by  its 
equivalent  y  from  x  =  g  *(y) .  Thus 
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p2(y)  *  p(x)  =  p (y) 


-A  <  y  <_  B 


-A 

B 

is 


Let 


-  ■  /.: 


p(x)dx  =  Prob  [x  <  -A] 


(U.8) 


/CO 

p(x)dx  =  Prob  [  X  >.  B] 


When  A  <  B,  clearly  >  a2  •  Also 


Prob[-A  <_  x  _<  B]  =  J  p(x)dx  =  1  -  -  a. 


(b.9) 


The  quantities  and  a 2  represent  for  the  output  record 


Thus 


*i  =  1  p2(y)dy  =  Prob  [y  i  -a] 

*2  -  ^  p2(y)dy  *  Prob  fv  >_  b] 

f  B 

1  *  al  "  a2  =  J  p2Cy)dy  =  Prob  [-A  <  y  <  B] 


(fc.10) 


(U.ll) 


Note  that  for  A  <  B,  one  has 


'p(-A)  >  p(B) 


(k.12) 


Consider  now  the  shape  of  p2(y).  As  y  approaches  the 
value  —A  from  inside4  [y  >  —A]  , 


P2(-A  )  =  p(-A) 


(fc.13) 


However,  as  y  approaches  the  value  —A  from  outside^  [y  <  -A], 
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d.U) 


P2f-A”J  - 

For  an  assumed  a  <  1,  there  results 

P2(-a")  >  P2(-A+)  (1.15) 

The  reverse  inequality  would  be  true  if  a  >  1. 

Similarly,  as  y  approaches  the  value  B  from  inside, 

[y  <  b]  , 

P2(B")  =  pCB)  (1.16) 

while  as  y  approaches  the  value  B  from  outside  [y  >  B] , 

p2(B  +  )  =  (1.17) 

Here,  for  an  assumed  b  >  1,  these  results 

P2(B+)  <  p2(B‘)  (1.18) 

with  a  reverse  inequality  if  b  <  1.  This  treatment  shows  that 
the  shape  of  P2(y)  under  the  assumption  that  0  <  9.  <  1^  b  >  !_> 
must  be  as  sketched  in  Figure  1.3.  For  the  assumed  values  here, 

observe  that  «  > 

1  2 

P2(-a”)  >  p2(B+)  (1.19) 


-105- 


Suppose  the  three-slope  system  is  hard-limited  at  x 


This  gives  the  output  probability  density  function 


P2  (y)  =  P  (y)  -A  _<  y  <  B 

=  a^5 (y  +  A)  y  =  -A 

=  a 2 <5  (y  -  B)  y  =  B 


where  6 ( y )  is  the  usual  delta  function.  Assuming  that 


Figure  ^.5  Output  PDF  for  hard-limited  three-slope  system 


- . 2  CATENARY  SYSTEMS 


A  second  example  of  a  nonlinear  nonsymmetrical  system  is  a 
catenary  system  as  sketched  in  Figure  1.6.  Let 


y  =  g (x)  =  -aA' 


x  <  -A 


=  a  ( x  +  2 Ax) 


-A  <.  x  <.  B 


=  a  (B  +  2AB) 


x  >  B 


Figure  K6  Catenary  system. 


For  definiteness,  assume  here  that  the  constants 


A  >  B  >  0, 


a  >  0 


[i 


The  inverse  relations  are 

x  =  g  1  (y)  =  any  value  of  x  <_  -A 

2 

when  y  =  -aA 

=  -A  +  Va2  +  (y/a)  f  -aA2  <  y  <  a(B2  +  2AB)  (it. 25) 

=  any  value  of  x  B 
2 

when  y  =  a(B  +  2 AB) 

The  middle  range  is  obtained  by  solving  the  quadratic  equation 

x2  +  2  Ax  —  (y/a)  =0  (1.2  6) 

where  only  the  positive  square  root  operation  makes  sense.  The 
derivative  relations  are 

=  g,(x)  =  0  x  i  "A 

=  2a  (x  +  A)  -A  <  x  <  B  (1.27) 

=  0  x  >  B 

Thus,  the  derivative  is  discontinuous  at  x  =  B. 

Let  p(x)  be  the  input  probability  density  function  and 
P2<y)  be  the  output  probability  density  function.  For  this 
problem,  from  Equation  1.5, 

p7(y)  =  -  (1.28) 

|dy/dx| 

only  in  the  range  where  -A  <  x  <  B,  and  where  the  value  of  x 
in  the  right-hand  side  should  be  replaced  by  its  equivalent  y 
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for  y  in  the  range 


P\~A  +  Va  +  (y/a) 


2a  VA  +  (y/a) 


This  has  the  form  of  a  chi-square  probability  density  function 
when  p(x)  is  assumed  to  be  Gaussian.  Different  relations  are 
required  for  x  <_-A  and  x  >_  B . 

Assume  the  input  probability  function  p(x)  is  Gaussian 
with  ux  =  0  and  =  1  as  sketched  in  Figure  b.7  where 


p(x)  =  — —  exp(-x2/2) 


TZZfl  l 


Figure  4.7  Gaussian  input  PDF 


The  assumed  values  for  A  and  B  are  such  that  A  >  B  >  0.  Let 


Clearly, 


f  -A 

a  =  I  p(x)dx  =  Prob  [x  <_  -A] 
J  —00 


00 

8=1  p(x)dx  =  Prob  [x  _>  B] 
B 


(Ml) 


a  <  B 


(M2) 


Also, 


rB 

Prob  [-A  <  x  <  B]  =  /  p(x)dx  =  1  —  a  - B 

J  -  A 


(M3) 


All  of  the  values  of  x  <_  -h  are  associated  with  the  single 

2 

value  of  y  equal  to  y  =  -aA  .  Hence  the  output  probability 
density  function 


P2  (y)  =  a  <5  (y  +  aA  ; 


y  =  -aA 


(MM 


where  6(y)  is  the  usual  delta  function.  This  gives 


Prob  [y  =  -aA^J  =  a 


(J*.35) 


Similarly,  all  the  values  of  x  >  B  are  associated  with  the 

2 

single  value  of  y  =  a(B  +  2AB) .  Hence  the  output  probability 
density  function 
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P2 (y)  =  8  6  (y  -  a[B2  +  2AB])  at  y  =  a(B2  +  2AB) 
where  6(y)  is  the  usual  delta  function.  This  gives 

Prob  [y  =  a(B2  +  2AB)]  =  8 

.  2 
There  is  zero  probability  here  that  y  <  -aA  or  that 

2 

y  >  a(B  +  2AB) .  The  probability  that  y  falls  inside  these 
two  bounds  is 

Prob  £-aA2  <  y  <  a  (A2  +  2AB)J  =  1  -  a  -  8 

A  typical  shape  for  p2(y)  is  sketched  in  Figure  l*.8. 

Observe  that  P2(y)  approaches  infinity  as  y  approaches  the 

i  2  2 

value  -aA  from  inside  [y  >  -aA  ] ,  but  p2(y)  stays  finite 
as  y  approaches  the  value  a(B2  +  2AB)  from  inside. 


Figure  U.8  Output  PDF  for  catenary  system. 
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Input/Output  Spectrum  Relations  for  Catenary  System 


In  the  range  from  -A  _<  x  <_  B ,  the  catenary  output  y(t) 
is  assumed  to  consist  of  two  components: 

(1)  a  linear  term  that  is  proportional  to  x(t) , 

2 

(2)  a  nonlinear  term  that  is  proportional  to  x  (t) . 

In  terms  of  constants  k^  and  k2  , 

y  ( t)  =  k1x(t)  +  k2x2(t)  (U. 39) 

This  equation  can  be  made  to  apply  to  more  general  situations 
as  shown  in  Figured. 9  by  letting  the  catenary  output  be  represent¬ 
ed  by 


y(t)  =  yx  ( t )  +  y2(t)  +  n(t) 


(U.  1*0) 


where 


y^(t)  =  linear  output  due  to  x(t) 

2 

y2(t)  =  nonlinear  output  due  to  v(t)  =  x  (t) 

n(t)  =  uncorrelated  zero  mean  Gaussian  output  noise 

In  Figure  1*.9  the  constant  parameter  linear  system  frequency  response 
functions  H^(f)  and  H2(f)  incorporate  the  constants  k-^  and 
k2,  respectively.  Note  that  Figure  U.9  is  a  special  case  of  Figure  1.8 
that  can  be  analyzed  by  formulas  in  Section  1.5*  Figure  1*.9  can  be  replaced 
by  a  tvo-input/single-output  linear  model  with  inputs  x(t)  and  v(t)  that  will 
be  uncorrelated  when  x(t)  is  Gaussian. 
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Figure  ^.9  Catenary  system  model. 

In  the  frequency  domain,  the  governing  relation  of  this 
catenary  system  model  is 

Y  ( f )  =  Y1(f)  +  Y  2  ( f )  +  N(f)  (U.  Ul) 

where  Y^(f),  Y2(f)  and  N(f)  are  Fourier  transforms  of  y^Ct), 


y2(t)  and 

n  (t)  , 

respectively. 

In 

terms  of  X(f)  and 

V(f)  , 

the  Fourier 

transforms  of  x(t) 

and 

v ( t )  =  x2 (t) , 

Y 

i(f>  = 

H1(f)  X(f) 

(h.h2) 

Y 

2(f)  = 

H2(f)V(f) 

V 

(f)  = 

fm  vi 

J  —  00 

dt  = 

[  X  (el)  X  (f  -  *)d* 

— oo 

(h.hh) 

Note  that  V(f)  can  be  computed  by  two  different  methods,  either 
directly  from  v(t)  or  indirectly  from  X(f) . 

For  Gaussian  input  data  x(t)  as  assumed  here,  the  output 
terms  y^(t)  and  y2(t)  will  be  uncorrelated.  Hence,  the  two- 
sided  output  spectral  density  function  S^y(f)  satisfies  the 
relation 
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syylf)  ‘  sy1y1'£)  +  sy,y,(f)  +  snn<f> 


y2y2 


(U.U51 


In  terms  of  Sxx(f)  and  Svv(f) ,  the  two-sided  spectral  density 

2 

functions  of  x(t)  and  v(t)  =  x  (t) , 


Vi*'5  -  iH1(f>isxx(f) 


(U.U6) 


sy0y  <f)  =  lH2(f>  I  svv<f> 


(U.UT) 


Svv(f)  =  T  E[V*(f)V(f)] 


(h.bQ) 


Also j  fo r  f  &  0, 


Svv(f)  =  2 


Sxx(«)Sxx(f  -  a>da 


(b.hg) 


Note  that  S  (f)  can  be  computed  by  two  different  methods, 
either  directly  from  V(f)  or  indirectly  from  S  (f)  . 

XX 

Given  H-^(f)  and  H2(f),  plus  measurement  only  of  x(t), 

the  above  equations  show  how  to  determine  the  spectral  properties 

of  y, (t)  and  y~(t),  namely,  S  (f)  and  S  (f).  If 

Ylyl  y2y2 

y(t)  is  measured  as  well  as  x(t) ,  then  one  can  also  determine 

the  two-sided  spectral  density  function  S  (f)  of  the  unmeasured 

nn 

n(t)  by  the  relation 


S _ (f)  =  s  (f)  -  S„  „  (f)  -  S„  ..  (f) 


ylyl 


y2y2 


(u.50) 


A  good  model  occurs  when  S  (f)  is  small  compared  to  S  (f) . 

nn  yy 

For  situations  where  H^(f)  and  H2(f)  are  not  known,  opti¬ 
mum  estimates  of  their  properties  can  be  obtained  that  minimize 


Snn(f)  provided  one  can  make  simultaneous  measurements  of  x(t) 
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and  y(t).  This  system  identification  problem  can  be  solved  using 
the  relations 


Sxy(£)  ”  sxyi(f)  “  Hl'f)sxx(£> 


(It. 51) 


Svy(£)  =  H2(£)Sw(f) 


(1.52) 


Thus 


H1(f) 


H2(f)  - 


S  (f ) 
xy  ' 

Sxx(f> 

Svy(f> 

Svv<f> 


(fc.53) 


(h.5b) 


Here,  S  (f)  is  the  two-sided  cross-spectral  density  function 
between  x(t)  and  y(t),  and  Svy(f)  is  the  two-sided  cross- 
spectral  density  function  betw*'  v(t)  and  y(t).  These 
quantities  can  be  computed  directly  by 


Sxy(f)  "  f  E[X*(f)Y(f)] 


(fc.55) 


Svy(f)  =  k  Etv*(f)Y(f)] 


(l*.56) 


The  quantity  S  (f)  can  also  be  computed  indirectly  by 


Svy(£> 


-r 

j  —  o 


Sxxy(c,'£  -  °>da 


(**.57) 


where  S  (a,f  -  a)  is  the  two-sided  cross-bispectrum  (also 

AXy 

called  the  second-order  cross-spectral  density  function)  defined 
for  all  f  and  g  by  the  two-dimensional  function 
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Sxxy(f'g)  =  f  E[X*(f)X*(g)Y(f  +  g)  ]  (fc.58) 

The  integrand 

Sxxy(ct'f  "  a)  =  f  E[X*(a)X*(f  -  o)Y(f)]  (14.59) 

is  a  special  case  of  sXxy^f,g^  when  f  =  ot  and  g  =  (f  -  a). 

For  f  ^  0,  one  can  also  compute  S  (f)  from  Equations  1.19 
and  b.52  by  the  relation 


r 

j  — o 


Svy(£)  =  2H2(f)  I  _  SXXMSXXU  -  Ode. 


(U.60) 


Hence,  for  f  ^  0,  the  two-sided  cross-bispectrum  satisfies 


Sxxy (ot,f  "  a)  =  2H2(*)Sxx(a)Sxx(f  "  a)  (U-6l) 

This  provides  an  alternate  way  to  determine  H2(f)  from  compu¬ 
tation  of  the  other  quantities.  Consider  a  change  in  variables 
where  a  =  g  and  f  =  2g.  Then  (f  -  a)  =  g  so  that  one  obtains 

Sxxy(g'g)  =  2H2(2g)Sxx(g)  {k-62) 

Now,  replacing  g  by  f  gives 

W£'£>  *  2H2‘2£>Sxx(£>  (1,-63> 

Solving  for  H 2 ( 2 f )  shows 

H  (2f)  =  f  x21y  (f-,f-  (h.6h) 

2S2  (f) 
xx 
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Finally,  replacing  f  by  (f/2)  proves 


H2(f) 


S  (f/2,f/2) 
2Sxx(f/2> 


(U.65) 


This  concludes  the  material  on  system  identification  and  on 
input/output  spectrum  relations  for  the  catenary  system. 
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4.3  SMOOTH- LIMITER  SYSTEMS 

Consider  the  nonlinear  nonsymmetric  smooth-limiter  system 
sketched  in  Figure  L.10  where 


y  =  g( 


exp(-t  /2)dt 


x  >  0 


(It. 66) 


=  0  x  <  0 

Equation  U.66  is  the  equation  for  the  normal  probability  integral. 


y 


Figure  1*.10  Smooth -limiter  system. 

For  this  nonlinear  transformation,  special  values  are 


g(0)  =  0 
g(l)  =  0.683 
g ( 2)  =  0.954 
g ( 3)  =  0.997 
g  (°°)  =  1.000 


(it. 67) 
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Tables  are  available  that  give  g(x)  as  a  function  of  x,  and 
that  can  also  be  used  to  find  x  for  any  g(x)  in  range 
0  <  g(x)  <  1. 

The  inverse  relations  are 

x  =  g  ^(y)  =  any  value  of  x  <_  0  for  y  =  0 

=  unique  value  of  x  in  range  0  <  x  <  »  for  0  <  y  <  1 

=  no  value  of  x  for  y  <  0  or  for  y  >  1 

The  derivative  relations  are 

=  g*(x)  =  y/Jexp(-x2/2)  x  >  0 

(U.68) 

=  0  x  <  0 

Here  the  derivative  is  discontinuous  at  x  =  0. 

Assume  that  p(x) ,  the  input  probability  density  function, 
is  Gaussian  with  mean  zero  and  standard  deviation  a  where 

p(x)  =  — - —  exp(-x2/2o2)  (U .  69) 

aV7\ r 

All  of  the  values  of  x  £  0  are  associated  with  the  single  value 

of  y  =  0.  Hence  P2(y),  the  output  probability  density  func¬ 
tion,  is  such  that 

P2(y)  -  \  6  (y)  at  y  =  0  (U.70) 

where  the  factor  (1/2)  occurs  because  the  probability  is  (1/2)  that  xi  0, 
There  are  no  values  of  x  where  y  <0  or  y  >  1.  Hence 

p2(y)  =  0  for  y  <  0  or  y  >  1  (U.71) 
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In  the  range  where  0  <  y  <  1 


from  Equation  1.5*  the  result  is 


p2(y>  -  p(x)  =  y,r — ^ — 1  0.-72) 

|  dy/dx  |  Z  |_exp(-xV2)  _ 

where  the  value  of  x  on  the  right-hand  side  should  be  replaced 
by  its  equivalent  y  from  x  =  g  1(y).  As  y  .approaches  the 
value  0  from  the  inside  where  y  >  0  ,  the  value  of  x  also 
approaches  zero  so  that 

p2«>+)  -  yf  [eti]  -  h  {u-m 

As  y  approaches  the  value  1  from  the  inside  where  y  <  1,  the  value  of 
x  approaches  infinity  so  that 

p2(0=  « 

A  typical  shape  for  p2(y)  is  sketched  in  Figure  U.u. 


(U.7M 


Figure  H.ll  Output  PDF  for  smooth -limiter  system. 

It  is  of  interest  to  be  able  to  determine  the  output  auto¬ 
correlation  function  and  the  input/output  cross-correlation  func¬ 
tion  for  Gaussian  input  data  through  this  smooth -limiter  system. 
Such  formulas  can  be  found  using  Price's  Theorem  of  Equation  1.6 
and  Bussgang's  Theorem  of  Equation  1.7. 
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Output  Autocorrelation  Fuhttion  for  Smooth-Limiter  System 


By  Price's  Theorem,  the  output  autocorrelation  function 
Ryy { t )  and  input  autocorrelation  Rxx(t)  f°r  Gaussian  input 
data  satisfies  the  equation 


3-V(T) 

9Rxx,T> 


=  E[g,(x1)g*  (x  )  ] 


s  |  |  g  (xL)g  (x2)p(x1,x2)dXidx2 


(U.75) 


where 


x1  =  x(t) 


X2  =  X ( t  +  x) 


y1  =  y(t) 
y2  =  y ( t  +  x) 


U.76) 


The  correlation  functions 


Rxx(T)  *  Etxlx2] 


Ryy(T)  -  E[ylY2) 


U.T7) 


Derivative  values  from  Equation  k.6 8  are 


g1(x1)  =  yf  exp (-x2/2)  for  x  >  0,  otherwise  zero 
g,(x2)  =  exp(-x2/2)  for  x  >  0,  otherwise  zero 


U.78) 


The  joint  probability  density  function 


p(xltx2)  = 


"22 

^  ”^xl  +  x2  “  2P^1x2) 

_______  exp  - 

v4  -  p2  (_  2°2(1  _  p2) 


( h. 79) 
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Thus 


3R  (x) 

yy 

3R  (T) 
XX 


irVvTT?  'o 


a 


exp 


-  (x2  +  x2  -  2px1x2)  (x2+x2) 


2  .  2, 


2oT(1  -  p2) 


dx^dj 


(U. 80) 


The  bracketed  quantity  coital  nr 


(xl  +  x2  -  2px1x2)  (x2  +  x2)  (x2+x2)[l+02(l-p2)]  -  2px^x. 


2o2 ( 1  -  p2) 


2  2 

2cj  ( 1  -  pz) 


x2  +  x2  -  2pXlx2/[l  +  c2(1  -  p2)] 
2a2  ( 1  -  p2)/[l  +  o2(l  -  p2)] 


2  2 

X1 +  x2  "  2poxlx2 


2o2  (1  -  p2)/[l  +  a2  (1  -  p2)  ] 


X1  +  x2  ~  2poXlx2 
2a2 (!  -  p2) 


(Ml) 


where  pQ  and  aQ  are  defined  by 


po  " 


[1  +  a2 ( 1  -  p2)] 


R  (x) 

With  p  =  p  (T)  =  - 


XX 


(M2) 


2  2 
°o(1-po> 


=  a2  (1  -  p2)/[l  +  o 2  ( 1  - 


P2)  ]  =  a2  (1  -  p2)  (po/p) 


a  =  a 


P0(l  -  p‘> 


P(  1  -  p‘> 


(k.83) 
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By  these  substitutions, 


?Ryy(T)  .  i  ["[ 

?Rxx<T>  J°J 


2  2 

-  ^  +  x2  -  2pQx1x2) 

2  2 
20^(1  -p2) 


dx1dx2  (U.S1*) 


The  integral  is  now  in  the  form  of  a  joint  Gaussian  probability 
density  function  over  one  quadrant  where 


2  2 

n,  -(x,  +  x2  -  2p  x1x2)  , 

-  exp  -  dx  dx0  =  T 

,  2  tf. - T  ,  2,,  2,  12  4 

2tto  VI  -  p  2a  (1  -  p  ) 

o  o  o  o 


"12  "  4 


(**.85) 


Hence 


3Rw(T)  .  (n!ZLiA 

3Rxx(T>  V2”  02vfT7 


p-  V'l  -  p‘ 

fi\— _ 

(2tt*  , - i 

p  vi  -  7: 


(**.86) 


p2  = 


[1  +  a2 ( 1  -  p2)]2 


( 1*  -  87 ) 


2  [1  +  o2(l  -  p2)  ]2  -  p2  1  +  2o2(l-p2)  +  o4(l-p2)2-p2 

[1  +  a2  (1  -  p2)]2  (p/pQ)2 


=  ( PQ/P) 2  (1  -  P2)  (1  +  o2)2  -  a4p2 


(I4.88) 
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This  gives 


° e I1  -  °2> 

p2d  -  »i) 


1 

r ,,  L  2.2  42, 

[  (1  +  o  )  -op] 


1 

[(1  4-  O2)2-  R*  (T)] 


( 1* .  89 ) 


Now 


3V(T) 


3Rxx<T>  2*[<1  +  a2)2  -  R2  (t)Jl/2 


Thus 


Ryy(T)  “fe)sin"1 


Rxx(t) 

rr^ 


+  c 


Qsing  the  fact  that  the  indefinite  integral 


/ 


dx  .  -1/  x  .  ,  „ 

=  sin  f— 1+  C 


J~2 - 7 

va  -  x 


(*)■ 


where  C  is  a  constant  of  integration.  At  t  =  0, 
Ryy (0)  =  ^y  =fe)sin  1 


*"  a2  “> 


l  +  o‘ 


+  C 


(U.90) 


(1*.  91) 


(U. 92) 


( U . 93 ) 


from  which  C  can  be  determined. 
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Input/Output  Cross-Correlation  Function  for  Smooth-Limiter  System 


By  Bussgang's  Theorem,  the  input/output  cross-correlation 

function  R  (t)  for  Gaussian  input  data  satisfies  the  equation 
xy 


Rxx(T)  f" 

R  ( T )  =  - 5 —  /  xg  (x)  p  (x)  dx 

xy  a  J0 


(k.9h. 


where 


g  (x)  =  |/|  J  exp  (-t2/ 2)  dt 


x  >  0 


(1.95: 


p(x)  =  — - —  exp(-  x2/2 cr2) 

a/2i 


(1.96) 


ntegration  by  parts  formula  gives 


dv  =  uv  -  /  v  du 


-I 


(1.97) 


u  =  g (x) 


dv  =  xp(x)dx  (1.98) 


Then 


du  =  g *  (x)  dx  = 


exp(-  x  /2)dx 


(1.99) 


V  =  -  o  p  (x) 


(1.100) 


At  x  =  0  and  at  x  =  00 ,  the  product  uv  =  0.  Hence 
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g  (x)  p  (x)  dx 


p(x)  g 


(x)  dx 


where  is  defined  by 


a,  = 


1  +  a' 


The  integrand  is  now  in  the  form  of  a  first-order  Gaussian 
probability  density  function  for  x  >  0  satisfying 


Hence 


and 


g {x; p (x) dx  = 


/2ir(l  +  a2) 


R . (t)  = 


Rxx(t> 


\^2v 


(1  +  az) 


This  concludes  Section  1. 


(1.101) 


(1.102) 


(1.103) 


(1.101) 


(1.105) 
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